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Podstawowe Oznaczenia 


Wszystkie wektory — jeżeli tego nie zaznaczono wyraźnie — są wektorami kolumno¬ 
wymi 


d 

H lub K 
P lub N 

M lub S 

x = (ai,... ,x d ) T 


x = (1, x±,.. .,x d ) T 


xW = (4 n) ,...,4 n) ) r 

x (n) = (l,x { { l) ,...,x < f ) ) T 
y = (yii • • • dm) t 


liczba cech, składowych wektora danych (u staty¬ 
styków oznaczana symbolem p, w MatLabie: R , 
u Osowskiego N) 
liczba neuronów warstwy ukrytej 

liczba wektorów uczących (liczebność próbki), od 
pattern (u statystyków: rozmiar próbki n lub N) 
liczba wyjść sieci neuronowej, czyli liczba neuronów 
ostatniej warstwy sieci, w MatLabie S 
wejściowy wektor danych (wektor kolumnowy) 
(u statystyków: osobnik, jest pisany jako wektor 
wierszowy), często pisze się górne wskaźniki bez na¬ 
wiasów okrągłych 

wejściowy wektor danych w którym jako pierwszy 
element podstawiono tożsamościowo jedynkę (wte¬ 
dy wektor wag otrzymuje jako pierwszą składową 
wartość wq oznaczającą stałą bo (bias) neuronu 

wejściowy wektor danych (wektor kolumnowy) 
o numerze n 

(u statystyków: n-ty osobnik, pisany jako wektor 
wierszowy); w MatLabie oznaczany jako wektor p 

wektor danych dla p-tego wzorca, powiększony 
o element „1” stojący na pierwszym miejscu 
wyjściowy wektor wyników, wyprodukowanych 
przez sieć na podstawie wektora wejściowego x. 
Pisząc bardziej dokładnie, mamy y = y(x). 
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y(n) = 




t = 


t” 


(t 


n 4.n \T 

1 > * * • l M) 


W i = {Wi i,w i2 ,---,w ij ) T 


(fc) 

u. ■ 

31 


W# x ( d+ i) 


to samo, co poprzednio, ale napisane explicite dla 
n-tego osobnika. 

Mamy oczywiście: y ^ = y(x^) 

docelowy wektor wyników (target). Sieć - na pod¬ 
stawie wektora wejściowego x ma wyprodukować 
wynik y = (y i,..., j/m) T możliwie zbliżony do war¬ 
tości wektora docelowego t = (ty,... £m) T ; w przy¬ 
padku sieci z jednym wyjściem wektor wartości do¬ 
celowych redukuje się do skalaru 
to samo, co poprzednio, ale napisane explicite dla 
n-tego osobnika, czyli dla n-tego wzorca 
waga z jaką i-ty neuron sumuje j- tą cechę, lub ak¬ 
tywację j-e go neuronu poprzedniej warstwy, i = 
j = 0,l,...,d 

wektor wag i-tego neuronu, przy założeniu, że i-ty 
neuron otrzymuje sygnały od j innych neuronów; 
w zależności od kontekstu może to być również 
wektor zawierający na pierwszym miejscu sta¬ 
łą Wio = ho, czyli bias i-tego neuronu: w; = 
(w l0 ,m ź i,tc,: 2 , • • • ,Wij) T 

oznacza wagę połączenia pomiędzy j - tym 
neuronem warstwy fc-tej oraz i-tym neuronem 
warstwy k — 1-szej 

tablica wag, w której wiersze i oznaczają numery 
neuronów w rozważanej warstwie (zakładamy tu, 
że warstwa ta zawiera H neuronów) 
kolumna zerowa zawiera tzw. wagi zerowe (obcią¬ 
żenia, bias) 

kolumny 1,..., d oznaczają numery cech w danych 
wejściowych 


w Hx(d+1) 



w 


T 

H J 


W Hxd, b// x i - tablica wag o znaczeniu podobnym, jak w poprzed¬ 
nim punkcie, jednak bez wag zerowych, które wy¬ 
odrębniono w postaci wektora kolumnowego b 
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Ui lub Vi - wynik sumowania bodźców dochodzących do neu¬ 
ronu nr i od neuronów (wejść) poprzedniej war¬ 
stwy; kombinacja liniowa obliczana przez i-ty neu¬ 
ron na podstawie wczytanego wektora danych x 
wynosi odpowiednio: 

d d 

Ui = WijXj + Wio lub Ui = u>ijXj + bi 
3=1 3 = 1 

W MatLabie wynik u,; jest oznaczany jako rij. 
Przyjmując tożsamościowo Xq = 1 oraz bi = w^, 
wynik Ui możemy zapisać prościej jako 

d 

Ui = WijXj = wfx 
3=0 


/(•) - Funkcja aktywacji. Może być: skokowa (f. Heavisi- 

de'a, inaczej hardlimit), liniowa, sigmoidalna (lo¬ 
gistyczna lub tangensoidalna), typu softmax. 




Rozdział 1 


Wstęp 


1.1. Plan zajęć 

1. Wprowadzenie w tematykę sieci neuronowych, trochę historii, neurony biolo¬ 
giczne i synapsy, określenie perceptronu. 'Neural Network Toolbox’ MatLab’a, 
Obliczenia z MatLabem. 4+4 g. 

2. Sieć jednokierunkowa jedno- i wielowarstwowa, funkcje aktywacji i znajdywanie 
wag, diagramy Hintona. 4+4 g. 

3. Pakiet Netlab. Przykłady realizacji zagadnień regresyjnych lub aproksymacyj¬ 
nych; porównanie z metodami statystycznymi. 4+4 g. 

4. Sieci neuronowe o radialnych funkcjach bazowych, wybór liczby funkcji bazo¬ 
wych, stabilność rozwiązania. 4+4 g. 

5. Przykłady zastosowań sieci neuronowych w klasyfikacji i diagnostyce. 4+4 g. 

6. Samoorganizujące się sieci Kohonena, odwzorowanie punktów z wielowymiaro¬ 
wych przestrzeni z zachowaniem topologii konfiguracji, wizualizacja Sammona, 
pakiet SOMĆPAK. 4+4 g. 

7. Przykłady zastosowań sieci Kohonena do wizualizacji przebiegu procesu pro¬ 
dukcyjnego, w marketingu usług turystycznych i ocenie stabilności banków. 
Inne typy sieci neuronowych i przykłady ich zastosowań. 


Pracownia będzie oparta na pakiecie MatLab; będą wykorzystywane: 
MatLab, Neural Network Toolbox, 

NETLAB (autorzy: łan Nabney i Ch. Bishop), 

http://www.ncrg.aston .ac.uk 

SDM Toolbox (autorzy: J. Vesanto i inn.) 

http://www.cis.hut.fi/projects/somtoolbox 


2 


Rozdział 1. Wstęp 


1.2. Linki do plików z danymi: 

http://www.ics.uci.edu/pub/machine-learning-databases/ 

http://www.phys.uni. torun.pl/kmk/ 

http://citeseer.nj.nec.com/prechelt94proben/ 

http://kdd.ics.uci.edu/databases/tic/tic.html 

O sieciach Kohonena: 

http://www.student.ii.uni.wroc.pl/~Vol/ 



Rozdział 2 


Od biologicznego neuronu do 
perceptronu i madaline 

2.1. Elementy neurobiologii. Komórka nerwowa i jej 
struktura 

Mózg i system nerwowy nie stanowią struktury ciągłej, ale składają się z około try¬ 
liona (10 18 ) komórek, z czego około 100 miliardów (10 11 ) stanowią komórki nerwowe 
połączone w sieci (Korbicz i inn. [']), dzięki którym realizowane są funkcje inteligen¬ 
cji, emocji, pamięci i zdolności twórczych. 

Ciało komórki nerwowej (inaczej: neuronu) jest dość podobne do komórek in¬ 
nych tkanek, wyróżnia się jednak wielkością otaczających ciało wypustek w postaci 
rozkrzewionych gałązek, tzw. dendrytów. 

Z komórki nerwowej wychodzi długie włókno, nazywane aksonem, które na ogół 
rozgałęzia się w postaci tzw. drzewka aksonowego. Akson, rozwidlając się, dociera do 
wielu komórek, niemniej sygnał wyjściowy jest identyczny dla wszystkich odbiorców. 
Zakończenia gałązek aksonu stykają się z dendrytami innych neuronów, a miejsce 
styku nazywa się synapsą. 

Podstawowe zadanie neuronu sprowadza się do przyjmowania (poprzez dendryty), 
przetwarzania i dalszego przekazywania (poprzez akson) informacji w postaci bodź¬ 
ców elektrycznych. W fizjologii pobudzanie aksonu określa się jako wszystko albo nic. 
Oznacza to, że dostatecznie silny bodziec powoduje każdorazowo tę samą reakcję, 
zbyt słaby bodziec nie wyzwala żadnej reakcji. Każdy nadchodzący synapsą bodziec 
dochodzi do ciała komórkowego. 

Przewodzenie poprzez synapsy następuje zawsze tylko w jednym kierunku. 
Informacja wzdłuż wypustek (aksonów, dendrytów) jest przenoszona w postaci im¬ 
pulsów elektrycznych, nazywanych potencjałami czynnościowymi. 
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Rozdział 2. Od biologicznego neuronu do perceptronu i madaline 


ab c 

Zakończenia 


przedsynaptyczne 



Rysunek 2.1: Schemat biologicznego neuronu na podstawie książki Michajlika i Ra- 
motowskiego [5], str. 372, i pracy A. Kotuli [9]. Synapsy neuronu to: 

A - synapsa aksonowo-dendrytowa, B - synapsa aksonowo-somatyczna, 

C - synapsa aksonowo-aksonowa bliższa (zazwyczaj hamująca), 

D - synapsa aksonowo-aksonowa dalsza (zawsze hamująca). 


Mózg 1 , stanowiący centrum ludzkiego organizmu, jest częścią mózgowia obej¬ 
mującą półkule mózgowe i część wzrokową podwzgórza. Mózgowie przeciętnie waży 
u człowieka około 1.3 kg. Mózg pokryty jest mocno pofałdowaną warstwą kory mó¬ 
zgowej. Warstwa ta ma grubość około 3 mm i powierzchnię 2500 cm 2 3 . Kora mózgowa 
składa się z upakowanych gęsto komórek nerwowych o różnej wielkości i kształcie, 
tworzących kilka warstw. Liczbę neuronów szacuje się na 10 miliardów. Przyjmują 
one i wysyłają impulsy o częstotliwości 1 — 100 Hz, czasie trwania 1 — 2 ms, napięciu 
IOOtos szybkości propagacji 1 — 100 — . Liczba połączeń między komórkami szacowana 
jest na 10 15 . Tadeusiewicz w książce [15], str. 13, podaje, że szybkość pracy mózgu 
oszacować można na 10 18 £PynL£21 Dj a porównania 2 , firma IBM planuje skonstruować 
w ciągu najbliższych pięciu lat komputer Blue Gene, który byłby w stanie wykonać 
10 15 operacji na sekundę. Komputer Blue Gene byłby w takim przypadku około dwa 
miliony razy potężniejszy niż dostępne dzisiaj komputery osobiste, tysiąc razy po¬ 
tężniejszy od komputera Deep Blue 1 i pięćset razy potężniejszy od najszybszych 
obecnie komputerów na świecie. 

ha podstawie informacji podanych przez A. Kotulę [9] 

2 Kotula [9] — na podstawie informacji zawartych na stronach internetowych firmy IBM, dostęp¬ 
nych pod adresem http: //www. ibm. com. 

3 Na podstawie stron internetowych firmy IBM, Deep Blue jest komputerem specjalnie zaprojek¬ 
towanym do gry w szachy przez grupę specjalistów (Feng-Hsiung Hsu, Murray Campbell, Joe Hoane, 
Jerry Brody oraz C.J. Tan) pracujących dla firmy IBM. Projektowanie Deep Blue rozpoczęto w roku 
1989, ale już od 1985 Hsu zajmował się tym problemem. W roku 1997 Deep Blue rozegrał pierwszy 
słynny mecz z najlepszym wówczas szachistą świata, Garrim Kasparowem, w meczu rewanżowym 
w maju 1997r wygrał Deep Blue — por. http://www.research.ibm.com/deepblue/. 









2.2. Matematyczny model neuronu wg McCullocha i Pittsa 
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2.2. Matematyczny model neuronu wg McCullocha 
i Pittsa 

McCulloch i Pitts w roku 1943 jako pierwsi zaproponowali znacznie uproszczony 
w stosunku do rzeczywistego model neuronu, który do dziś jest podstawą większości 
modeli. Schemat takiego „matematycznego” działania neuronu jest podany na ry¬ 
sunku 2.2. 



Rysunek 2.2: Schemat działania neuronu o numerze i — według McCullocha i Pittsa 

Na rysunku tym mamy zaznaczony jeden neuron — ma on umownie numer i. Do 
neuronu tego zbiegają się sygnały (bodźce) - jest ich j. Neuron je sumuje — z wa¬ 
gami wn, Wi 2 ,. ■., Wij odpowiednio. Gdy obliczona wartość sumy przekroczy pew¬ 
ną wartość progową wm, specyficzną dla danego neuronu, następuje jego „zapłon”, 
inaczej mówiąc, neuron ten znajdzie się w stanie pobudzenia. Matematycznie stan 
pobudzenia neuronu wyraża się dwiema wartościami: 0, gdy pobudzenie neuronu nie 
przekroczyło jego specyficznej wartości progowej, i 1, gdy jest przeciwnie. 
Spróbujmy teraz te fakty zapisać matematycznie. Będziemy rozpatrywać neuron o nu¬ 
merze i ze specyficzną wartością progową w^. Założymy, że stan pobudzenia neu¬ 
ronu jest zjawiskiem dyskretnym zmieniającym się w czasie r w stałych odstępach 
czasu Ar. Oznaczmy wartość pobudzenia i-go neuronu w czasie r symbolem ą(t). 
Oczywiście wartość neuronu w chwili r + Ar zależy od tego, jak były pobudzone (do¬ 
starczające mu bodźce) neurony z jego otoczenia oznaczone tu umownie jako zbiór 
{i/} — w momencie poprzedzającym moment r. Neuron oblicza sumę ważoną dostar¬ 
czanych mu sygnałów. O ile po dodaniu do wyznaczonej przez neuron sumy wartości 
progowej otrzyma się liczbę dodatnią, następuje zapłon. Wyrazić to można w sposób 
następujący: 



( 2 . 1 ) 


j 


Zmienna -i(t) może mieć wartość 1, gdy i -ty neuron znajduje się w chwili r w stanie 
zapłonu, lub 0, gdy tak nie jest. 

Wagi występujące w powyższym wzorze odzwierciedlają istotność synapsy łą¬ 
czącej neuron ż-ty i j- ty. Wagi mogą przyjmować zarówno dodatnie jak i ujemne 
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Rozdział 2. Od biologicznego neuronu do perceptronu i madaline 


wartości: 

Wij 


> 0 
= 0 
< 0 


odpowiednik synapsy pobudzającej 

odpowiednik braku połączenia pomiędzy neuronami 

odpowiednik synapsy hamującej. 


Natomiast funkcja ©(a) występująca we wzorze McCuloclra i Pittsa to funkcja He- 
aviside’a ( hardlimit ) określona następująco: 


0(a) 


1 dla a > 0 
0 dla a < 0 


McCulloch i Pitts wykazali, że przy odpowiednio dobranych wagach Wj : j synchronicz¬ 
ny zespól takich neuronów może wykonać te same obliczenia, co uniwersalna maszyna 
licząca. 


Dalsze prace poszły w kierunku: 

• użycia innych funkcji aktywacji - umożliwia to modelowanie procesów nieli¬ 
niowych, 

• przedstawienia sygnału z t nie jako procesu dyskretnego, ale jako procesu cią¬ 
głego. 

Stosuje się następujące uogólnienie modelu McCullocha i Pittsa: 

Zi = g(^2 W ij z j + w io)- (2.2) 

3 

We wzorze tym nie uzależnia się stanu pobudzenia neuronu od czasu r. Funkcja 
progowa ©(•) zastąpiona jest przez funkcję g (-), zwaną funkcją aktywacji (funkcją 
wygładzającą, funkcją przejścia, funkcją wzmocnienia). Model (2.2) uwzględnia aktu¬ 
alizację Zi w dowolnej chwili, umożliwia nieliniowość, Zi jest funkcją ciągła określającą 
stan neuronu w chwili t. 


2.3. Perceptrony 

Co to jest perceptron 

Pojęcie perceptronu 4 zostało wprowadzone przez Rosenblatta 5 i przez długi czas było 
bardzo popularne, przede wszystkim w zagadnieniach związanych z rozpoznawaniem 
wzorców, lub — używając języka innego środowiska — zagadnień dyskryminacji i kla¬ 
syfikacji danych. 

4 u Osowskiego perceptron nie jest specjalnie omawiany, chociaż na str. 22 mówi się o uczeniu 
bezgradientowym na zasadzie reguły perceptronu', pojęcie perceptronu jest omawiane specjalnie np. 
w książce Hertza i wsp., Rozdział 5, lub Korbicza i wsp., Rozdział 2 

5 Rosenblatt F., 'Principles of Neurodynamics’, Washington D.C., Spartan Press 1962, również 
ten sam autor, 'The perceptron: A probabilistic model for information storage and organisation in 
the brain’, Psychological Review, Vol. 65, 1958, str. 386-408 
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Percepton Architecture 

Input Percepton Layer 


V_ ) 



a = hardlim( W X + b ) 


The Perceptron Leaming Rule 

W”" = W°" + eP 


where 


b ""' 1 = b" u + e 

e = t - a 


Rysunek 2.3: Perceptron z R wejściami i S wyjściami 


Rozpatrujemy sieć neuronową o następującej architekturze: R wejść, S neuronów 
w warstwie ukrytej, funkcja aktywacji typu hard limit (Heaviside'a), liczba wyjść 
równa liczbie neuronów warstwy ukrytej (=K). Na wyjściu mogą pojawiać się tylko 
zera lub jedynki. Schemat takiej sieci jest przedstawiony na rysunku 2.3. 

Funkcja aktywacji typu hardlimit jest przedstawiona na rysunku 2.4. 

Idea perceptronu jest zawarta w następujących zasadach (por. Korbicz i wsp. [8]): 

1. Perceptron jest siecią zbudowaną ze sztucznych neuronów; każdy sztuczny neu¬ 
ron - którego budowa jest przedstawiona na rysunku 2.5 — jest aktywowany 
funkcją hardlimit typu unipolarnego (zwracającą wartości 1 lub 0) lub bipolar¬ 
nego (zwracającą tylko wartości 1 lub —1). 

2. Sieć perceptronową można podzielić jednoznacznie na ściśle uporządkowane 
i rozłączne klasy elementów, zwanych warstwami, wśród których wyróżnić moż¬ 
na warstwę wejściową i wyjściową. Pozostałe noszą nazwę warstw ukrytych. 
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Y 

+1 


-1 

a = hardlim(n) 

Hard Limit Transfer Function 


Rysunek 2.4: Funkcja aktywacji typu hardlimit , unipolarna 


Perceptron Neuron 


Input 


PerceptronNeuron 



Where... 

R = # Elements 
in input vector 


Rysunek 2.5: Perceptron z R wejściami, 1 neuronem warstwy ukrytej i 1 wyjściem 


3. Perceptron nie zawiera połączeń między elementami należącymi do tej samej 
warstwy. 

4. Połączenia pomiędzy warstwami są asymetryczne i skierowane zgodnie z ich 
uporządkowaniem, tzn. od warstwy wejściowej do pierwszej warstwy ukrytej, 
następnie od pierwszej do drugiej warstwy ukrytej, itd. aż do warstwy wyjścio¬ 
wej. Nie ma połączeń zwrotnych. 

Konsekwencją tych założeń jest, że każdy neuron — niezależnie od pozostałych 

uczy się rozpoznawać pewną strukturę. 

Popatrzmy na neuron o numerze s, (s = 1,..., S). Neuron ten najpierw oblicza 
sumę u s = 6^+^^ w S iXi , która następnie zostaje przetransformowana przez funkcję 
aktywacji typu hardlimit (por. rys. 2.4). Wynik transformacji, czyli wartość g(u s ), 
zostaje skierowana bezpośrednio na wyjście. 

Stała b s jest tzw. zerowym obciążeniem ( bias ). 
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2.4. Perceptron prosty — Klasyfikacja do dwóch klas 

Nie ma dodatkowego przetwarzania przez neuron warstwy wyjściowej, czyli na 
wyjściu pojawiają się bezpośrednio wartości g(u a ) — przyjmujące wartości 1 lub 0, 
odpowiednio. 

Tak zdefiniowana sieć — perceptron nadaje się znakomicie do zagadnień klasyfi¬ 
kacyjnych, gdy naszym celem jest zaklasyfikowanie wektora danych [x\,... ,a;#] do 
jednej z S klas danych. Działanie perceptronu można w tym przypadku interpretować 
następująco: 

Każdy neuron „specjalizuje” się w rozpoznawaniu jednej, np. s-tej klasy danych. 

W trakcie uczenia neuron przystosowuje swoje wagi tak, aby móc rozpoznać prób¬ 
ki (tj. wektory danych) należące do s-tej klasy. „Obce” próbki (tzn. nie należące do 
s-tej klasy) są przez ten neuron odrzucane. 

Wtedy, na zakończenie procesu uczenia, wektor wagowy w s = [u> s i,..., ui s _r] wraz 
ze stałą b s tworzy kombinację liniową zmiennych wejściowych xi ,..., xr czyli wartość 
■u s = b s + ZZi w si x i =b s + wjx, taką że; 

• duże wartości u s świadczą za tym, że przedstawiany wektor x jest podobny do 
elementów s-tej klasy; 

• małe wartości u s świadczą za tym, że przedstawiany wektor x jest niepodobny 
do elementów s-tej klasy; 

(Własności te wynikają z własności funkcji hardlimit ) 

Perceptron tak zdefiniowany potrafi nauczyć się rozpoznawać próbki należące 
do różnych klas pod warunkiem, że klasy te są liniowo separabilne, czyli 
rozłączne. 

Proces uczenia („reguła perceptronu”) jest bardzo prosta. Każdy neuron (s) po¬ 
winien znaleźć taki wektor wag w s i stalą b s , które by pozwoliły na efektywne rozpo¬ 
znanie struktury danych należącej do klasy s. Proces znajdywania odpowiednich wag 
odbywa się bezgradientową metodą iteracyjną omówioną w następnym podrozdziale. 

Rozenblatt pokazał, że jeśli klasy są liniowo separabilne, to zaproponowany przez 
niego algorytm potrafi znaleźć wektory wagowe pozwalające na poprawne rozpozna¬ 
nie elementów należących do poszczególnych klas. 

2.4. Perceptron prosty — Klasyfikacja do dwóch klas 

W przypadku dwóch klas wystarczy jeden neuron warstwy ukrytej, który będzie 
się specjalizował w rozpoznawaniu próbek pochodzących z klasy pierwszej. Schemat 
takiej sieci jest pokazany na rysunku 2.5. Oblicza się wtedy jedną wartość u = bo + 
W\X\ + ■ ■ ■ + wrXr. Duże wartości u świadczą za przynależnością próbki do klasy 
pierwszej, a małe wartości za przynależnością do klasy drugiej. Granicą rozdzielającą 
obie klasy jest hiperpłaszczyzna L 


L : b 0 + w\x i + • • • + wrXr = 0. 


(2.3) 
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Rozdział 2. Od biologicznego neuronu do perceptronu i madaline 


Zauważmy, że wektor w tworzący hiperpłaszczyznę (2.3) jest wektorem ortogonal¬ 
nym do tej hiperpłaszczyzny leżącym w półprzestrzeni w której powinna znajdować 
się rozpoznawana klasa. 

W przypadku tylko dwóch zmiennych wejściowych x\ i X 2 (czyli: R = 2) hiper- 
płaszczyzna (2.3) redukuje się do prostej bo + W\X\ + W 2 X 2 = 0. Prosta ta rozdziela 
płaszczyznę ( 21 , 2 : 2 ) na dwie części: tam, gdzie y > 0 oraz gdzie y <= 0. Wartości 
dodatnie świadczą o przynależności próbki do klasy 1, A wartości zerowe lub ujemne 
wskazują na przynależność do klasy 2. Sytuacja taka jest pokazana na rysunku 2.6. 

Decision Boundary 


L 



Rysunek 2.6: Perceptron określony dla 2 zmiennych x\ i 22 . Granica decyzyjna (de¬ 
cision boundary) dla określania przynależności do jednej z dwóch klas danych jest 
określona prostą L o równaniu w x+6 = 0, gdzie W = [wpi 11 : 1 , 2 ]. Prosta ta rozdziela 
elementy klasy 1 i klasy 2 


2.5. Uczenie (trenowanie) perceptronu — reguła per¬ 
ceptronu 

Zasadę tę wyrażono w opisie MatLab NN Toolbox za pomocą frazy, że uczenie sieci 
jest dokonywane „by altering the vector w towarcl input vectors to be classified as 1 
and away from vectors to be classified as 0”. Zasada ta jest zilustrowana na rysunku 
2.7. 












2.5. Uczenie (trenowanie) perceptronu — reguła perceptronu 
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Uczenie każdego neuronu przebiega niezależnie od pozostałych. 

Sieć uczymy (trenujemy) na podstawie tzw. próbki uczącej. Celem jest wyznacze¬ 
nie wag, na podstawie których można, dla danego wektora x, wyznaczyć wartości y 
możliwie zgodne z podanymi wartościami docelowymi t (Przypomnijmy, że percep¬ 
tron wytwarza na wyjściu jako y tylko wartości 1 lub 0; natomiast wartości docelowe 
(t) podawane w próbce uczącej, powinny określać, czy chcemy otrzymać na danym 
wyjściu wartość 1 czy też 0, co można deklarować conajmniej na 2 sposoby). 

Przypuśćmy, że sieć ma tylko jedno wyjście, czyli M = 1. Niech e oznacza błąd 
sieci: 

e = d — y. 

Zauważmy, że błąd ten może przyjmować tylko trzy wartości: 0, 1 i -1. 

Zasada uczenia (trenowania) sieci jest następująca: 

Przyjmujemy (np. losowo) jakiś wektor w jako wektor wag, i jakąś początkową 
wartość obciążenia b. 

Przedstawiamy sieci (nowy, lub też powtórnie, stary) wektor x dla którego jest 
znana wartość docelowa t. Obliczamy błąd e = t— e. W zależności od wielkości błędu 
obliczamy poprawkę Aw i Ab rozważając trzy przypadki: 


Xi 



x { 



Rysunek 2.7: Poprawianie wag w trakcie uczenia. Lewa: Wektor x nie został wykazany 
jako należący do grupy 1 (d — y = e = +1). Wobec tego wektor w należy dosunąć 
do wektora x obliczając nowy wektor wag jako \v new = w oid + x Prawa: Wektor x 
został rozpoznany błędnie jako należący do grupy 1 (d — y = e = —1). Wobec tego 
wektor wag w należy odsunąć od wektora x odejmując ten ostatni od wektora w. 
Uwaga: Fakt , że wyrysowany wektor x jest U do w - jest całkowicie przypadkowy. 
W rzeczywistości kąt między tymi dwoma wektorami może być zupełnie dowolny 


Przypadek 1. Obliczona wartość błędu jest równa zero, czyli e = t — y = 0. Oznacza 
to, że prognoza y była dobra, a wiec nic nie musimy zmieniać. Praktycznie 
mamy tu: Aw = 0 i Ab = 0. 
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Przypadek 2. Obliczona wartość błędu jest równa 1, czyli e = t — y = 1. Oznacza to, 
że sieć prognozowała y błędnie. Powinna była zaprognozować „1” a zaprognozo- 
wała „0”. Można stąd wnioskować, że wagi nie są dostatecznie zaadaptowane do 
tego, żeby rozpoznawać wektory klasy „1”. Wobec tego robimy taką poprawkę, 
żeby zbliżyć wektor wag w do wektora x. 

Praktycznie: dodajemy wektor x do wektora w, co można zapisać w postaci 
reguł: Aw = x i Ab = 1. 

Przypadek 3. Obliczona wartość błędu jest równa — 1: e = t — y = — 1. Oznacza to, 
że sieć prognozowała y błędnie. Powinna była prognozować „0” a prognozowała 
„1”. Można stąd wnioskować, że wagi zbyt często prognozują „1”, nawet gdy 
mamy do czynienia z wektorami klasy przeciwnej. W tej sytuacji robimy po¬ 
prawkę mającą na celu odsunięcie wektor wag w od aktualnie przedstawianego 
sieci wektora x. 

Praktycznie: odejmujemy wektor x od wektora w, co można zapisać w postaci 
reguł: Aw = x i A6 = — 1. 


Nasze postępowanie jest zilustrowane na rysunku 2.7, gdzie pokazano dodawanie 
i odejmowanie wektora x do (od) wektora w. Wymienione trzy przypadki można 
zapisać razem w postaci jednego wyrażenia: 

Aw = (t — y)x T = ex T , oraz Ab = (t — a)(l) = e, 


lub 

w new = w oW + ex, b new = b old + e. (2.4) 

We wzorze powyższym symbol w oznacza wektor wagowy dowolnego neuronu. 
Zarówno w jak i x są wektorami kolumnowymi. 

Jeżeli mamy sieć składającą się z S neuronów w warstwie ukrytej, to wszystkie 
wektory wagowe możemy zestawić w macierz W o wymiarach S x R: 


W S xR 



(2.5) 


Wtedy regułę uczenia całego perceptronu (czyli wszystkich jego wektorów Wago¬ 
wych) możemy zapisać za pomocą jednego równania: 

W new = W oid + ex T ) 
h new = h old + e 

gdzie W = W Sxfi , e = (ei,.. .,e s ) T , no i oczywiście x = (xi, .. .,x R ) T . 








2.6. Kiedy model perceptronowy daje dobre wyniki 
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2.6. Kiedy model perceptronowy daje dobre wyniki 

Przedstawiona metoda (dyskryminacji liniowej) jest względnie prosta i jest dość sze¬ 
roko używana w pakietach komercyjnych przeznaczonych do analizy danych pomia¬ 
rowych. 

Należy jednak zaznaczyć, że metoda ta działa dobrze, gdy klasy są rozłączne, 
czyli mają własność liniowej separabilności linear separability. 

Wtedy rzeczywiście wystarczy prosta lub hiperpłaszczyzna żeby rozdzielić po¬ 
szczególne klasy, a sieć neuronowa jest zdolna nauczyć się w skończonej liczbie kro¬ 
ków skonstruować hiperpłaszczyznę rozdzielającą (czyli wyznaczyć wektor wag, na 
podstawie którego taką hiperpłaszczyznę można określić). 

Zasada perceptronu została dość zdecydowanie i pesymistycznie zaopiniowana 
w książce M. Minsky i S. Pappert, Terceptrons’, Cambridge Ma, MIT Press, 1969. 
Pokazali oni, że perceptrony Rosenblatta mają ograniczoną zdolność uczenia się, 
a niektórych struktur nie potrafią się wręcz nauczyć. Opinia ta wpłynęła na zahamo¬ 
wanie — na okres conaj mniej 10 lat — badań naukowych związanych ze sztucznymi 
sieciami neuronowymi. 

Jednak jest znanych wiele przypadków, gdy analizowane dane nie posiadają wła¬ 
sności „liniowej separabilności”, i nie są to bynajmniej dane wydumane, ale wystę¬ 
pujące w praktyce. 

W przypadku takich danych metoda perceptronową, tak jak ją przedstawiono 
wyżej, nie będzie dawać dobrych wyników i należy sięgnąć po inne metod, np. sko¬ 
rzystać ze sztucznej sieci neuronowej o wzbogaconej strukturze (np. sieci zawierające 
więcej neuronów lub wielowarstwowe). Innymi, bardziej efektywnymi metodami mogą 
również okazać się metody statystyczne, np.: 

• Liniowa funkcja dyskryminacyjne Fishera 

• Kanoniczne funkcje dyskryminacyjne 

• Postępowanie oparte na ilorazie wiarogodności 

• Modelowanie nieparametryczne (kernel density) 

• Modele skomponowane z mieszanin rozkładów (mixture models) 


2.7. Moduły MatLaba demonstrujące działanie per¬ 
ceptronu 

MatLabowski Pakiet NN Toolbox oferuje następujące moduły demonstracyjne obra¬ 
zujące działanie perceptronu: 


1. nnd4db, Decision Boundaries 

2. nnd4pr, Perceptron learning rule. Pick boundaries 
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3. demopl, Classification with a 2-input perceptron 

4. demop4, Outlier input vector 

5. demop5, Normalized perceptron rule 

6. demop6, Linearly non-separable vectors. 

Wszystkie moduły ilustrują działanie klasycznego perceptronu dla przypadku, 
gdy liczba cech (wymiar wektora wejściowego) wynosi 2. (W MatLab-ie wektor wej¬ 
ściowy jest oznaczany jako p, a jego składowe jako p(l) i p(2) — zamiast naszych 
oznaczeń X\ i X2', wartości docelowe są oznaczane jako t - od ang. słowa target ). 

W przypadku dwu-wymiarowych danych wejściowych granicę decyzyjną (czyli 
prostą y = WiX\ + 1x2X2 + b) można zilustrować na płaszczyźnie. 

Moduły nnd4db i nnd4pr demonstrują uczenie się perceptronu - który w zależno¬ 
ści od podanych danych wyznacza odpowiednie wagi Wi,W2 oraz obciążenie ( bias ) b. 
Moduły te pozwalają na interakcyjne określanie danych wejściowych oraz manualne 
korygowanie prostej decyzyjnej. Moduły te pochodzą z książki Hagana i wsp. [3]. 

Na uwagę zasługuje moduł demop4: pokazuje on, jak mocno odstająca obserwacja 
może utrudnić uczenie się sieci. 


2.8. Adaline i Madaline 

Mniej więcej w tym samym czasie co Rosenblatt (to jest około r. 1960), Widrow 
wraz z współpracownikami zbudował dość podobną sieć co perceptron, ale opartą 
na liniowej funkcji adaptacji, która dopuszczała na wyjściu dowolne liczby rzeczywi¬ 
ste. Zbudowana przez Widrowa sztuczna sieć składała się z elementów nazywanych 
Adaline, czyli Adaptive Linear Element. Cala sieć jako wielokrotność (multiple) tych 
liniowych elementów „adaline” została nazwana Madaline. 

Do tak zbudowanej sieci Widrow z współpracownikami określili regułę uczenia 
opartą na minimalizacji sumy kwadratów błędów, czyli metodzie najmniejszych kwa¬ 
dratów. 


2.9. Typologia sieci ze względu na pewne jej własno¬ 
ści formalne 

2.9.1. Podstawowe architektury sieci neuronowych 

Rozróżniamy następujące typy sieci neuronowych 6 


6 Wiadomości głównie na podstawie książki Osowskiego [12], str. 18 
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1. Sieci jednokierunkowe jednowarstwowe i wielowarstwowe. 

Warstwa wejściowa nie liczy się. 

Sieci jednowarstwowe używane są w praktyce bardzo rzadko (przykłady: per¬ 
ceptron prosty, adaline). 

Najczęściej używa się sieci dwuwarstwowych. Sieci takie składają się z war¬ 
stwy wejściowej, warstwy ukrytej o H neuronach, i warstwy wyjściowej o M 
neuronach. 

Przykłady sieci jednokierunkowych jedno- i dwuwarstwowych są pokazane na 
rysunku 2.8. 



Rysunek 2.8: Przykłady sieci jednokierunkowych jedno- i dwuwarstwowych 


Sieć dwuwarstwowa potrafi rozpoznać klasy obiektów, które geometrycznie są 
ulokowane w obszarach będących kombinacją iloczynów i sum mnogościowych 
stanów aktywnych. Przykładowo, dla perceptronu stan aktywny neuronu wy¬ 
znacza półprzestrzeń; dla kilku neuronów te półprzestrzenie przecinają się i two¬ 
rzą zbiory wielościenne wypukłe (por. Korbicz i inn. [8], str. 41, Tadeusiewicz 
[16] str. 107-108). 

Sieci o więcej niż 2 warstwach używa się przy modelowaniu skomplikowanych 
i wybitnie nieliniowych związków między danymi wejściowymi i wartościami 
docelowymi (targetem). 

2. Sieci rekurencyjne, mogą być jedno- i wielowarstwowe, posiadają np. zdolność 
uwzględniania elementów opóźnienia przy analizie sygnałów. Znajdują szerokie 
zastosowanie w teorii sterowania. 

3. Sieci komórkowe'. Neurony znajdują się w specyficznej konfiguracji geometrycz- 


'Korbicz i inn. [8], str. 30; termin ten nie występuje u Osowskiego 
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Rozdział 2. Od biologicznego neuronu do perceptronu i madaline 


A 


sygnał 
wejściowy 



sygnał 

sprzężenia 

zwrotnego 


Rysunek 2.9: Przykład sprzężenia zwrotnego zachodzącego dla sygnałów przetwarza¬ 
nych przez jeden neuron — na postawie książki Tadeusiewicza [16], str. 244 

nej umożliwiającej określenie sąsiedztwa. Występują sprzężenia zwrotne między 
elementami przetwarzającymi, ale tylko między elementami znajdującymi się 
w „sąsiedztwie”. 

2.9.2. Metody uczenia sieci 

Wyróżniamy następujące metody uczenia: 8 

1. Uczenie pod nadzorem (z nauczycielem). 

Sieć najpierw — na podstawie próbki uczącej — uczy się rozpoznawać pewne 
klasy lub relacje. Sieć „uczy” się w ten sposób, że w miarę prezentacji wzor¬ 
ców jej wektory wagowe zostają zaadaptowane (zmodyfikowane) z myślą, aby 
zminimalizować błąd predykcji, czyli różnicę między wartością y prognozowaną 
przez sieć, a odpowiednią wartością docelową t podaną w próbce uczącej. 

Próbki (wzorce) przedstawiane sieci przy uczeniu mają postać 



gdzie 

wskaźnik n indeksuje (numeruje) wzorce próbki uczącej, 

liczba N jest liczebnością próbki uczącej, 

x" = [xi, x %,..., x]]] T - jest n-tym wektorem danych, 

t n = [i™,..., t'1 4 ] - jest wektorem wartości docelowych, odpowiadających ra-tej 
próbce, 

M - oznacza liczbę wyjść sieci. 

8 Osowski [ L2], str. 20-32, ze zmianą w oznaczaniu numeracji wzorców: u nas n, u Osowskiego p; 
i wartości docelowych: u nas t , u Osowskiego d 
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2. Uczenie z krytykiem. Jest to odmiana uczenia pod nadzorem. Sieć nie otrzymu¬ 
je informacji o wartościach docelowych t, ale jedynie informacje, czy podjęta 
przez system akcja (np. zmiana wartości wag) daje wyniki pozytywne czy też 
negatywne w sensie pożądanego zachowania systemu. 

3. Uczenie samoorganizujące się typu Hebba. Na zasadzie analogii do obserwacji 
neurobiologicznych (doświadczenia Pawłowa z bodźcami) przyjmuje się, że: 

• przy jednoczesnym stanie pobudzenia obu neuronów waga powiązań mię¬ 
dzy tymi neuronami wzrasta, 

• jeżeli dwa neurony nie są jednocześnie pobudzone, to waga połączeń mię¬ 
dzy nimi maleje. 

Uczenie jest procesem iteracyjnym, zmiana wag w kolejnych cyklach uczących 
może być zapisana następująco (k oznacza numer iteracji lub cyklu uczącego): 

Wij(k + 1 ) = Wij(k) + A Wij(k), gdzie A Wij(k) = F(xj,yi) 

przy czym F(xj 1 yi) jest funkcją stanu sygnału wejściowego Xj , zwanego presy- 
naptycznym, oraz stanu wyjściowego iji, zwanego postsynaptycznym. 

W klasycznym ujęciu Hebba jako funkcję F przyjmowano funkcję iloczynową 
obu sygnałów, a zmiana wag wyrażała się jako: 

A Wij(k) = r]Xj(k)yi(k), 

przy czym 77 oznacza tzw. współczynnik uczenia. 

Uczenie neuronu z zastosowaniem uczenia Hebba może odbywać się w trybie 
bez nauczyciela lub z nauczycielem. Znajduje zastosowanie przy nadzwyczaj 
ważnym zagadnieniu redukcji wymiarowości (kompresji) danych. Istotny wkład 
przy opracowywaniu algorytmów określających tzw. czynniki główne (składo¬ 
we główne Hotellinga, transformata Karhunen-Loeve’go, ang. Principal compo- 
nents) wnieśli T. Sejnowski, E. Oja i T.D. Sanger. 

4. Uczenie samoorganizujące się typu konkurencyjnego. Najbardziej znanym ty¬ 
pem sieci wykorzystującej ten sposób uczenia jest tzw. sieć Kohonena (SOM, 
Self Organizing Map). 

Specyficzną charakterystyką tego typu uczenia jest, że przy prezentacji nowe¬ 
go wzorca tylko jeden neuron może zostać pobudzony, a pozostałe pozostają 
w stanie spoczynku. Z tego powodu ten typ uczenia jest również określany jako 
WTA ( Winner Takes AU). Uczenie przebiega w ten sposób, że grupa neuronów 
współzawodniczących otrzymuje te same sygnały wejściowe x". i oblicza swoją 
aktywację Ui = wfx n . W wyniku porównania aktywacji m (wskaźnik i indek¬ 
suje tu neurony współzawodniczące) zwycięża ten neuron, którego wartość u,; 
jest największa. Neuron zwycięzca przyjmuje na swoim wyjściu stan 1 (i może 
np. zaadaptować swoje wagi), natomiast pozostałe (przegrywające) przyjmują 
stan 0 , a ich wagi pozostają niezmienione. 
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Rozdział 2. Od biologicznego neuronu do perceptronu i madaline 


2.9.3. Zdolności do uogólnienia 

Chcielibyśmy, aby sztuczne sieci neuronowe rozpoznawały również obiekty, czy też 
klasy obiektów, których nie pokazaliśmy w próbce uczącej. Niektóre sieci (np. sieć 
Fahlmana oparta o zasadę korelacji kaskadowej) potrafią rozpoznać, że przedstawiany 
wzorzec nie należy do żadnej dotychczas przez tę sieć rozpoznanych klas. 

Ocieramy się tu o zagadnienie „inteligencji” sztucznej sieci neuronowej. Do nie¬ 
dawna uważano, że zdolność uogólniania jest właściwa tylko umysłowi ludzkiemu. 



Rozdział 3 


Sieci neuronowe 

jednokierunkowe 

dwuwarstwowe 

3.1. Określanie architektury sieci 

Architekturę sieci jednokierunkowych określają następujące parametry: liczba neu¬ 
ronów warstwy wejściowej, liczba warstw ukrytych oraz liczba neuronów w każdej 
z tych warstw, liczba neuronów warstwy wyjściowej, funkcje aktywacji dla każdej 
warstwy ukrytej i warstwy wyjściowej. 

Generalnie sieć może wykonywać dwie funkcje: 

1 . uczyć się, tj. adaptować swoje wagi przy zadanej swojej architekturze 0, 

2 . wykonywać swoje właściwe zadanie, tj. na podstawie swojego algorytmu wyko¬ 
nywać predykcję, czyli obliczać wartości y(x; W, 0) 

Czynność 1 czyli uczenie sieci może być wykonywane wieloma sposobami przy użyciu 
różnych algorytmów. 

Czynność 2 jest już jednoznaczna. 

Przykład sieci jednokierunkowej dwuwarstwowej jest podany na rysunku 3.1. 

Architekturę tej sieci określają parametry © = (d, id, M ,/(.), (/>(.)) o następują¬ 
cym znaczeniu: 

d — liczba wejść, na rysunku: d = 2, 

H — liczba neuronów warstwy ukrytej, na rysunku: H = 3, 
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Rozdział 3. Sieci neuronowe jednokierunkowe dwuwarstwowe 


1 



v = f(u) 


Rysunek 3.1: Sieć jednokierunkowa dwuwarstwowa. Neuron i {i = 1,, H) pierwszej 
warstwy oblicza sumę Ui = = ^2y=i w ij x j> a następnie swoją aktywację Vi = 

podobnie neuron i z drugiej warstwy (i = 1 oblicza ważoną sumę 

dochodzących do niego bodźców: u[ 2 ' > = > a naa tępnie swoją aktywację 

Ui = <pK- 2) ) 


M — liczba neuronów warstwy wyjściowej, na rysunku: M = 2, 

/(•) oraz ip(-) — funkcje aktywacji warstwy ukrytej i warstwy wyjściowej 

Neurony warstw ukrytych i warstwy wyjściowej spełniają rolę sumatorów dochodzą¬ 
cych do nich - z warstw poprzednich — bodźców. Sumowanie jest ważone, w związ¬ 
ku z tym każdy neuron ma przypisany swój odrębny wektor wagowy, z jakim sumuje 
dochodzące do niego bodźce. 

W przypadku sieci dwuwarstwowej mamy dwa układy wektorów wagowych: 

wi, ..., w h ~ wektory wag przy neuronach warstwy ukrytej; w /, = [who, Whi, • • •, Whd] 
h = l,.... H (wskaźnik h oznacza numer neuronu danej warstwy; drugi wskaź¬ 
nik numer neuronu poprzedniej warstwy), 

w) ,..., - wektory wag przy neuronach warstwy wyjściowej; w j = [■ WjQ,Wj\ ,. .. , 

j = 1,... ,M, 

Wektory te można umieścić w macierzach i z wyraźnym zaznaczeniem 



3.2. Funkcjonowanie sieci po wytrenowaniu jej 
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— za pomocą górnego wskaźnika — której warstwy dane wagi dotyczą: 


w*- 1 ) — 

vv Hx(d+ 1) — 

' (wf , ) T ' 
(w<'>y 

— 

> vv M x (H+l) 

' (w‘ 2 y ■ 
(wfy 


. rfy. 


. (w2>r. 


Aby sieć funkcjonowała prawidłowo, tzn. wykonywała swoje zasadnicze czynno¬ 
ści (mogą nimi być: predykcja, klasyfikacja, grupowanie na zasadzie podobieństwa 
lub niepodobieństwa, filtrowanie i odszumianie sygnałów), powinny być określone 
wektory wag występujące w algorytmie wykonywanym przez sieć. Wagi te na ogól 
są wynikiem uczenia (trenowania) sieci, (w rozdziale 2 omówiliśmy już trenowanie 
perceptronu, dzisiaj omówimy jeszcze trenowanie prostej sieci zawierającej elementy 
adaline). 

3.2. Funkcjonowanie sieci po wytrenowaniu jej 

Celem trenowania jest otrzymanie takich współczynników wagowych (wektorów wag), 
aby po przedstawieniu sieci jakiejś (opcjonalnie nowej) próbki czyli wektora x sieć 
mogła dokonać predykcji lub też rozpoznania, do jakiej kategorii wzorców dana prób¬ 
ka x należy. 

Aby to zrobić, sieć wykonuje - dla każdego jej prezentowanego wektora x 
następujące obliczenia przebiegające w dwóch etapach: 

1. Najpierw każdy neuron warstwy ukrytej działa niezależnie jako sumator obli¬ 
czając ważoną sumę sygnałów wejściowych: 

Uh = w £x. 

Z obliczonej sumy Uh następuje wyznaczenie aktywacji neuronu Vh■ Do wyzna¬ 
czenia tej aktywacji używa się funkcji / zadanej przy określaniu architektury 
danej sieci: 

Vh = f(u h ), h=l,...,H. 

Obliczone w punkcie 1 aktywacje v \,..., vh są przekazywane do następnej war¬ 
stwy, w tym przypadku warstwy wyjściowej. 

2. Przekazywane przez poprzednią warstwę wartości Vh stanowią dane dla drugiej 
warstwy (wyjściowej), której neurony działają podobnie jak neurony warstwy 
ukrytej. A mianowicie: 

Każdy neuron j, (j = 1,..., M) wykonuje najpierw funkcję sumatora, oblicza¬ 
jąc sumę 

uf ] = (w] 2) ) T v, gdzie v = [«i,..., v H \- 
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Rozdział 3. Sieci neuronowe jednokierunkowe dwuwarstwowe 


Obliczony wynik sumowania jest następnie przetwarzany przez funkcję aktywa¬ 
cji ip(.) dostarczając ostatecznego wyniku yj pojawiającego się na j-tym wyjściu 
sieci: 

Vj = t P{ u< f)- 

3.3. Uczenie adaline — reguła delta Widrowa-Hoffa 

Adaline (Adaptive linear element) jest bardzo prostą siecią, podobną do percep¬ 
tronu prostego. Różni się funkcją aktywacji (jest nią funkcja signum lub funkcja 
purelin, czyli czysto liniowa) i sposobem uczenia. Zasady funkcjonowania adaline zo¬ 
stały opracowane przez Widrowa i współpracowników w latach 60-tych. Proste sieci 
adaline połączone w sieć jednokierunkową tworzą sieć Madaline, tj. Multiple adaline 
i sieci tego typu były używane do niedawna, a nawet podobno są używane jeszcze 
dzisiaj. 

Omówimy teraz bardziej szczegółowo zasady uczenia sieci typu adaline , czyli 
sieci z jednym wyjściem. 

Wprowadźmy następujące oznaczenia: 

x" = (x n i ,..., x n d) T - wektor danych (sygnał wejściowy), 

t n , y n - znana wartość docelowa (target) i wartość wyprodukowana przez sieć 
jako odpowiedź na sygnał x”. 

Różnica 

e n = t n — y n stanowi tzw. błąd sieci. 

e 2 = ( t n — y n ) 2 = (■ t n — y(x n ; w)) 2 stanowi kwadrat błędu. 

Chcielibyśmy, aby ten kwadrat błędu był możliwie mały. 

Jeżeli prezentacji było N (n = 1,..., AT), to chcielibyśmy wyznaczyć wektor wa- 
gowy w w ten sposób, żeby 

N 

'^2,e 2 n = mini 

n—1 

Czyli faktycznie stosujemy metodę najmniejszych kwadratów (least sąuares error). 
W zasadzie stosuje się dwie metody wyznaczania wag w: 

• Na zasadzie wyznaczenia wartości oczekiwanej kwadratu błędu E[e 2 ]. Wielkość 
błędu jest uważana za zmienną losową, a wartość oczekiwana E[e 2 ] jest wyzna¬ 
czana ze względu na rozkład prawdopodobieństwa zmiennych losowych (x, i). 

Jednak, aby obliczyć wartość oczekiwaną jakiejś zmiennej losowej, musimy znać 
rozkład prawdopodobieństwa tej zmiennej losowej. 

Metoda ta prowadzi do układu równań normalnych, których współczynnikami 
są kowariancje składowych wektora (x ±,..., Xd, t). 

1 Korbicz i inn., str. 47 
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Jeśli nie znamy rozkładu prawdopodobieństwa, a uważamy że posiadany zbiór 
wzorców jest reprezentatywny dla badanej populacji, to możemy przybliżyć 
wartość oczekiwaną sumy kwadratów błędów zaobserwowaną uśrednioną sumą 
kwadratów błędów, i szukać dla niej minimum. Również w tym przypadku 
dochodzimy do układu równań normalnych. 

Bardziej dokładną metodą — w przypadku znajomości rozkładu prawdopodo¬ 
bieństwa błędu — jest metoda największej wiarygodności. W przypadku gdy 
rozkład błędu jest rozkładem normalnym, otrzymujemy ten sam układ równań 
normalnych. 

• Na zasadzie iteracyjnego przybliżania wag. 

Prezentujemy sieci wiele razy (w różnej kolejności) wektory danych (nazywane 
również próbkami lub wzorcami). 

W każdym kroku k wagi ulegają adaptacji 

w (k + 1) = ~w(k) + Aw, 

przy czym poprawka Aw może być wykonywana na wiele różnych sposobów. 

Możemy również prezentować sieci cały zbiór wzorców i szukać minimum scho¬ 
dząc coraz niżej — aż do minimum. 

Jedna iteracja (jeden krok iteracyjny) jest nazywany epoką. 


W przypadku trenowania adaline Widrow i Hoff zastosowali następujące zasady: 
Określili błąd z pominięciem aktywacji, i tak określony błąd oznaczyli symbo¬ 
lem 5. Tym samym, kwadrat błędu czyli funkcja S 2 przy prezentacji n-tego wzorca 
jest określona jako: 

% = (*" - w T x") 2 . (3.1) 

Uaktualnianie wektora wag w w £:-tej iteracji, w której zaprezentowano sieci n-ty 
wzorzec x ra ,f" odbywa się wg już przedstawionej wyżej zasady 

w (fc + 1) = w (fc) + Aw, (3-2) 

przy czym jako A w przyjmuje się najczęściej kierunek najszybszego spadku funkcji 
błędu, a to oznacza kierunek przeciwny do gradientu funkcji błędu: 

Aw = -??V<5 2 (fc). (3.3) 

Współczynnik p w powyższym równaniu oznacza tzw. współczynnik uczenia - zaraz 
powiemy o nim coś więcej, najpierw jednak zajmiemy się gradientem funkcji błędu. 
Przy kwadracie błędu określonym wzorem (3.1) mamy 

V<5 2 — ~^—5 2 — —2 5 x" 

VO n — — zo n x- • 
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Rozdział 3. Sieci neuronowe jednokierunkowe dwuwarstwowe 


Stąd wynika wzór na poprawianie wektora w po kolejnej, k -tej iteracji: 


w (k + 1) = w (k) + 2p8 n x n . 


Reguła powyższa nosi nazwę reguły Widrowa-Hoffa; czasami jest nazywana 

regułą delty lub regułą najmniejszego błędu średnio kwadratowego. 

Wyznaczenie współczynnika 77 jest zadaniem skomplikowanym i wymaga znajo¬ 
mości statystyki prezentowanych wzorców wejściowych. Mówiąc ogólnie, współczyn¬ 
nik ten powinien być z zakresu: 


0 < 77 < 1/A 

max ? 

gdzie A max jest największą wartością własną macierzy A = A{xx T }. Jeżeli wzor¬ 
ce x są wystandaryzowane na (0,1), to macierz A jest macierzą korelacji między 
składowymi (cechami) danych wejściowych. 

Ogólnie, jeśli 77 jest zbyt duże, to proces poszukiwania optymalnego wektora jest 
rozbieżny; z kolei, jeśli wartość 77 jest zbyt mała, to zbieżność jest powolna. 

Nowsze algorytmy uzależniają wartość współczynnika 77 od czasu uczenia (numeru 
iteracji k) — w miarę upływu czasu uczenia 77 maleje. 

Istnieją dwie modyfikacje reguły Widrowa-Hoffa 2 

Pierwsza polega na tym, że nie dokonuje się modyfikacji wag po każdej prezentacji 
wzorca, ale dopiero po prezentacji pewnej ich liczby, np. Al. Po każdej pre¬ 
zentacji jednego wzorca x" oblicza się tylko poprawkę S n x n ; poprawki te się 
sumuje, a po prezentacji określonej liczby Al wzorców uśrednia się je, a na¬ 
stępnie — na podstawie takiej uśrednionej poprawki uaktualnia się wektor 
wag. W rezultacie otrzymuje się dokładniejszą aproksymację gradientu błędu 
średniokwadratowego Vó 2 . 

Druga modyfikacja polega na stosowaniu wzoru: 

w {k + 1 ) = w(fc) + 77(1 — p)( 5„x” + p(w(fc) — w (k — 1 )). 

Wprowadzony tu parametr p (p £ (0,1)) pozwala na zaniedbywanie zbyt du¬ 
żych fluktuacji wektora wag i spełnia rolę tzw. momentum znanego w teorii 
optymalizacji. 


2 Korbicz i inn. [8], str. 48 
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3.3.1. Odfiltrowywanie szumów za pomocą sieci zawierającej 
blok Tapped Delay Linę 

Na bazie struktury Adaline są zbudowane sieci wykonujące tzw. adaptacyjne filtro¬ 
wanie ( Adaptwe Filtering). 

W szczególności buduje się sieci zawierające tzw. Tapped Delay Linę. Na wejściu 
jest wprowadzany sygnał jednowymiarowy x(t) który jest zapamiętywany dla kilku 
(np. R) ostatnich momentów czasowych. Neuron sumujący zbiera R ostatnich sygna¬ 
łów, oblicza z nich ważoną kombinację liniową (wagami są wartości wag w i,..., wr), 
sumę przetwarza liniowo dodając do nich obciążenie wq , a wynik przekazuje na wyj¬ 
ście: 

R 

y(t) = purelin(w x + b) = Wix(t — i + 1) + wq. 

i=l 


3.3.2. Uczenie Madaline 

Sieć Madaline składa się z pojedynczych bloków adaline połączonych równolegle. 
Najczęściej stosujemy tu uczenie iteracyjne. Prezentowana sieci próbka ma postać 
(wskaźnik n oznacza tu indeks prezentowanej próbki, a wektor t oznacza prawidłowe 
wartości, które sieć powinna aproksymować): 

x”,t", gdzie x” = (x„i, ... ,x nd ) T , t" = (t nl , ..., t nM ) T ■ 

W rezultacie działania sieci otrzymujemy na M wyjściach M wartości 

Vi > ■ ■ • > 2/m i 

które zostały obliczone przez sieć z wczytanego wektora x” i aktualnych wag zapa¬ 
miętanych w wektorze w. 

Ponieważ wiemy, jakie powinny być wyniki, potrafimy określić wielkości błędów: 

lub ich kwadratów: 

i 6 ?) 2 = [*7 - y ?} 2 = [*? - y ( x n ; Wi )] 2 . 

We wzorach powyższych zaznaczono explicite, że odpowiedź sieci na j-tyrn wyjściu, 
czyli wartość y” jest funkcją wczytanej wartości x n oraz aktualnego wektora wag w j 
związanego z j-tym neuronem warstwy wyjściowej. 

Zarówno przy adaline jak i przy madaline korzystamy z reguły delta , która zakłada 
liniowe elementy przetwarzające, czyli że y” został wyznaczony ze wzoru: 

Vj =y(x";w i ) =wjx”. 
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Oznacza to, że wielkości delta-kwadrat są obliczane po prostu ze wzoru: 



gdzie j oznacza numer wyjścia. 

Niech E (od ang. Error ) oznacza ogólnie funkcję błędów, lub funkcję celu którą 
chcemy minimizować. Oznaczając błąd powstały na j-tym wyjściu symbolem Ej 
i sumując te błędy po wszystkich wyjściach otrzymujemy: 


M 


M 


M 



(3-4) 


Oczywiście E = E( w 1; ..., w M ) = Y.f=i E j( w j)- 


Z określenia błędu według wzoru (3.4) wynika, że zagadnienie minimalizacji wiel¬ 
kości E jako funkcji zmiennych wi,..., w m rozpada się na szereg niezależnych zadań 
minimalizacji poszczególnych wektorów w ;j (j = 1 ,..., M) charakterystycznych dla 
poszczególnych bloków j składających się z pojedynczych elementów adaline. 

Każdy z wektorów w j występuje tylko w jednym błędzie ó 7 , stąd zagadnienie 
sprowadza się do niezależnego obliczenia poprawki Awj dla każdego wektora w 
tak, jakby to był pojedynczy adaline. 

3.4. Jak określać funkcję błędu 

Widroff i Hoff przyjmowali, że należy minimizować sumę kwadratów błędów — wg 
zasady minimum sumy kwadratów. 

Reguła ta przetrwała do czasów dzisiejszych i funkcja błędu określona według tej 
zasady jest najbardziej popularnym kryterium również w czasach dzisiejszych. 
Jednak metoda ta w ostatnim dwudziestoleciu była poważnie krytykowana. 

Jako alternatywy proponowano 3 : 

• sumę wartości bezwzględnych błędów (norma L\) 

• bardziej specjalne funkcje, jak np. funkcja Hubera lub Hampela 

• mieszane funkcje, np. funkcja Karayiannisa 


E = l -\J2(tj - yj) 2 + (1 - A) - Vj)- 


M M 


3 por. Osowski [12], str. 41-44 
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W końcowej fazie ma być zapewniona minimalizacja wartości absolutnej błędu. 
Jako funkcję <f>i przyjmuje się 

$i(x) = h\[cosh((5x)\ 

która przy dużych wartościach j3 dąży do |x|. 


3.5. Moduły demonstrujące działanie sieci typu „Ada¬ 
line” 

(Tekst w głównej mierze oparty na dokumentacji MatLab Neural Networks version 3, Rozdział 4} 

Bardzo proste moduły demonstracyjne MatLaba pokazują, jak w najprostszym 
przypadku może komplikować się zagadnienie optymalizacji wektora wag. 

Moduły demolinl i demolin2 pokazują trenowanie sieci typu Adaline, tj. z liniową 
funkcją aktywacji. Rozpatruje się sieć która zawiera tylko jeden neuron. W rozpatry¬ 
wany przykładzie cala próba ucząca składa się z dwóch wzorców (osobników): p=[ 1.0 
— 1.2]; t=[ 0.5 1.0]; dla każdego osobnika mamy określoną tylko jedną cechę. 

Dla rozważanego problemu jest wykreślana powierzchnia błędu (error surface), 
oraz kontury tej powierzchni. 

Najlepszymi parametrami w i b są te, które wyznaczają minimum powierzchni 
błędu. Moduł demolinl znajduje te wartości za pomocą funkcji SOLVELIN, a moduł 
demolin2 metodą gradientów. 

Moduł demolin4 ilustruje działanie sieci w przypadku nieliniowego związku między 
zmienną wejściową x a zmienną docelową d. Przykładowo: 

P=[+1.0 +1.5 +3.0 -1.2]; 

T=[+0.5 +1.1 +3.0 -1.0]; 

W przypadku nieliniowej zależności błąd nigdy nie osiągnie wartości 0; jednak 
możemy zbliżyć się w rozsądnych granicach do przybliżonego rozwiązania. 

Dla wymienionego problemu ustala się liczbę iteracji (czyli maksymalną liczbę 
epok) net.trainParam.epochs= 15. 

Użytkownik podaje za pomocą funkcji ginput wartości początkowe un°), b Na¬ 
stępnie, posługując się reguła Widrowa i Hoffa, sieć „uczy się”, czyli wyznacza war¬ 
tości w, b dające coraz to mniejszy błąd dopasowania. Można zauważyć, że po kilku 
epokach błąd ten zaczyna się stabilizować. 

Moduł demolin5 ilustruje działanie sieci w przypadku, gdy układ parametrów nie 
jest określony w sposób jednoznaczny (ang. underdetermined problem), tzn. istnieje 
nieskończona liczba rozwiązań dla danego problemu. 

Przykład: p= [+1.0]; t= [+0.5]; Każda para wartości w,b spełniająca równość: 


b + wp = 0.5 
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jest rozwiązaniem dającym minimum błędu. Powierzchnia błędu (przedstawiona jako 
funkcja parametrów w i b) zawiera całą dolinę (rynnę) optymalnych rozwiązań. 

Moduł demolinó ilustruje działanie sieci w przypadku gdy przedstawiane sieci 
wzorce (próbki) są liniowo zależne, ale wartości docelowe nie odzwierciedlają tej za¬ 
leżności. 

Przykład: Próbki są dwuelementowe (x = {x\, £ 2 )), 1 mamy 3 takie próbki; zmien¬ 
na docelowa jest jedna: 

p=[ 1.0 2.0 3.0; 4.0 5.0 6.0]; t=[ 0.5 1.0 -1.0]; 

Tutaj w przedstawionym zestawie uczącym elementy środkowe (druga kolumna) 
w tablicy p są średnimi arytmetycznymi elementów kolumn 1 i 3; jednak zależność 
ta nie zachodzi elementu środkowego tablicy t. 

Znowu: nie można w tym przypadku zejść z błędem do zera; niemniej jednak sieć 
typu adaline — wykorzystująca przy uczeniu regułę Widrowa i Hoffa — znajduje 
przybliżone rozwiązanie. 

Moduł demolin7 ilustruje działanie sieci w przypadku uczenia się ze zbyt dużym 
współczynnikiem uczenia (Too large a learning ratę). 

Pokazane to jest na przykładzie: p=[+1.0 —1.2]; t=[+0.5 +1.0]; 

Dalsze dwa moduły ilustrują adaline sprzężoną z blokiem tapped delay linę. 

Moduł demolin8 (Adaptine linear layer ) pokazuje, jak sieć typu adaline z 1 jednym 
neuronem ukrytym i blokiem opóźniającym o 1 epokę potrafi nauczyć się nieliniowego 
związku między sygnałem wejściowym P(t) i sygnałem wyjściowym Y(t). W rezulta¬ 
cie, po pewnym czasie, predykcje sygnału wyjściowego y(t) pokrywają się z faktyczną 
wartością sygnału, a błąd predykcji jest praktycznie równy zeru. 

Moduł nndlOnc ilustruje działanie tzw. adaptacyjnego filtrowania (adaptine 
noise cancellation). Sieć składająca się z jednego adaline połączonego z blokiem 
opóźniającym o R epok potrafi z nakładających się dwóch sygnałów — właściwe¬ 
go i losowych zakłóceń, czyli tzw. białego szumu — odfiltrować zakłócenia. 

Moduł ten pokazuje również, jak szybko sieć uczy się prawidłowej predykcji, jeżeli 
w algorytmie uczenia zostaną uwzględnione: współczynnik uczenia i momentum. 

Sieć o podobnej architekturze może służyć do filtrowania (oczyszczenia z obcych 
sygnałów) np. transmisji głosu pilota przekazywanego w obecności szumu wytwarza¬ 
nego przez maszynę. Inny przykład, to odfiltrowywania echa (pogłosu) powstającego 
na łączach linii telefonicznych. 

Moduł nndlOlc powtarza demonstrację zaplanowaną przez Widrowa i Hoffa (por. 
[3], str. 10-37). Demonstracja ta jest oparta na 1 elemencie adaline, który potrafi 
nauczyć się rozpoznawać 3 litery (np. A, G, F). Próbka ucząca składa się z 3 liter, 
danych w postaci właściwej i przesuniętej. Kwadrat zawierający każdą literę został 
podzielony na 16 mniejszych kwadracików, z których niektóre zostały wyróżnione 
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kolorem zielonym. Tym samym każda próbka litery może być odczytana jako ciąg 
zero-jedynkowy, w którym „1” odpowiada pokolorowaniu na zielono. 

Przy danych wzorcach literek należy najpierw wytrenować sieć, wciskając przycisk 
train. Po wytrenowaniu można obejrzeć wagi wciskając przycisk weights - ukazuje się 
wtedy diagram Hintona. 

Jeżeli sieć wykształciła już swoje wagi (tzn. nauczyła się rozpoznawać 3 klasy 
znaków), to można kazać sieci rozpoznać literę podaną w polu „Test input”. Umiesz¬ 
czona na liczniku Widrowa-Hoffa wskazówka pokaże, do której z 3 wyuczonych klas 
znaków litera umieszczona w polu „Test input” jest najbardziej podobna. 



Rysunek 3.2: Demo nnlOlc, Learning classińcation lcl 
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— Niektóre instrukcje z~DEMHINT (Netlab) 

» nin=3; nhidden=5; nout=l; 

» net = mlp(nin, nhidden, nout, ’linear’); 

» net 
net = 

type: 'mlp* 
nin: 3 
nhidden: 5 
nout: 1 
nwts: 26 
actfn: ’linear’ 

wl: [3x5 double] 

bl: [-0.1532 -0.0432 -0.6800 -0.0486 -0.6285] 
w2: [5x1 double] 
b2: 0.2423 

» net.wl 
ans = 


-1.5304 

-0.3875 

-0.2240 

-0.5037 

-0.7243 

0.3973 

-0.0186 

-0.0054 

0.6721 

0.9509 

-0.0418 

-0.6402 

-0.2771 

-0.3750 

0.3447 


» net.w2 
ans = 

-0.2037 

0.3184 

0.3751 

0.0023 

-0.3220 

» [hl, h2] = mlphint(net); 

» hl 
hl = 1 

» h2 
h2 = 2 

» disp(’The MLP has been created with’) 

The MLP has been created with 
» disp([ J 1 int2str(nin) ’ inputs’]) 

3 inputs 

» disp([’ ’ int2str(nhidden) ’ hidden units’]) 

5 hidden units 

» disp([’ ’ int2str(nout) ’ outputs’]) 

1 outputs 
» disp(’ ’) 

» disp(’0ne figurę is produced for each layer of weights.’) 

» disp(’For each layer the fan-in weights are arranged in rows for each unit.’) 
» disp(’The bias weight is separated from the rest by a red vertical linę.’) 

» disp(’The area of each box is proportional to the weight value: positive’) 

» disp(’values are white, and negative are black.’) 
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3.5.1. Diagramy Hintona — proc. DEMHINT z Netlaba 

Pola kwadratów są proporcjonalne do wielkości wag. 

Czarne kwadraty — wagi ujemne; białe kwadraty — wagi dodatnie. 


neuron 1 


neuron 2 


neuron 3 


neuron 4 


neuron5 



x 3 




Bias 

b (1) 




Wagi: Wejście -- warstwa ukryta 


Rysunek 3.3: Wykres Hintona przedstawiający wagi pierwszej warstwy: nin=3, nhid- 
den=5 


nl n2 n3 n4 n5 b (2) 



Wagi: warstwa ukryta -- wyjście 


Rysunek 3.4: Wykres Hintona przedstawiający wagi dla drugiej warstwy : nhidden=5 
nout=l 


function demhint(nin, nhidden, nout) 

°/«DEMHINT Demonstration of Hinton diagram for 2-layer feed-forward network. 

°/o 

“/* Description 

# /o 

“/* DEMHINT plots a~Hinton diagram for a~2-layer feedforward network with 
°/* 5 inputs, 4 hidden units and 3 outputs. The weight vector is chosen 
1 from a~Gaussian distribution as described under MLP. 

1 

°/ e DEMHINT (NIN, NHIDDEN, NOUT) allows the user to specify the number of 
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7. inputs, hidden units and outputs. 

7. 

°/« See also 

7. HINTON, HINTMAT, MLP, MLPPAK, MLPUNPAK 

7 . 


7, Copyright (c) łan T Nabney (1996-9) 

if nargin < 1 nin = 5; end 
if nargin < 2 nhidden = 7; end 
if nargin < 3 nout = 3; end 

7, Fix the seed for reproducible results 

randn( 5 State *, 42); 

clc 

disp( 5 This demonstration illustrates the plotting of Hinton diagrams 5 ) 
disp( 5 for Multi-Layer Perceptron networks. 5 ) 
disp( 5 5 ) 

disp(’Press any key to continue. 5 ) 
pause 

net = mlp(nin, nhidden, nout, 5 linear 5 ); 

[hl, h2] = mlphint(net); 
clc 

disp( 5 The MLP has been created with 5 ) 
disp([ 5 1 int2str(nin) } inputs 5 ]) 

disp([ 5 5 int2str(nhidden) 5 hidden units 5 ]) 
disp([ 5 5 int2str(nout) 5 outputs 5 ]) 

disp( 5 5 ) 

disp( 5 0ne figurę is produced for each layer of weights. 5 ) 

disp( 5 For each layer the fan-in weights are arranged in rows for each unit. 5 ) 
disp( 5 The bias weight is separated from the rest by a~red vertical linę. 5 ) 
disp( 5 The area of each box is proportional to the weight value: positive 5 ) 
disp( 5 values are white, and negative are black. 5 ) 
disp( 5 5 ) 

disp( 5 Press any key to exit. 5 ); 
pause; 
delete(hl); 
delete(h2); 
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3.5.2. Przykłady grafiki z modułów demonstracyjnych: 



ff 

fc-- 








WeightW 


Override the randomly generated weight and bias by getting a mouse click on the 
Error Contour Plot wlth GINPUT: 

>> subplot(1,2,2); [net.lW{1,1)net.b(l)) = ginput(l X 

To show the path oł the training we will train only one epoch at a time and cali 

PLOTEP every epoch (codę not shown here). 

>> net = train(net,P,T>; 

The plot shows a history of the training. Each dot represents an epoch and the 
blue lines show each change madę by the learning rule (Wldrow-Hoff by default). 




Lewa: Demolin5, underdetermined problem. Powierzchnia błędu zawiera dolinę. 
Czerwone (puste, odizolowane) kółko oznacza rozwiązanie LSE. Plik lin5.eps. 
Prawa: Demolin7, Too large a learning ratę. Rozwiązanie odbija się od powierzchni 
błędu, zamiast schodzić do minimum. Proces szukania minimum przy zbyt dużym 
współczynniku uczenia może być rozbieżny. Plik lin7 



Demolin8, Sieć zawierająca na wejściu blok opóźniający sygnał potrafi nauczyć się 
prognozować nadchodzący sygnał z błędem równym zero. Plik lin8c 















































































Rozdział 4 


Zagadnienia dyskryminacji 
i klasyfikacji 


Przez zagadnienia dyskryminacji rozumiemy zagadnienie budowy funkcji dyskrymi¬ 
nacyjnych pozwalających na różnicowanie między kilkoma grupami danych. 

Zagadnienia klasyfikacji polegają na określaniu przynależności wektorów danych 
do poszczególnych klas, które są z góry zadane. 


4.1. Dyskryminacja i klasyfikacja dla dwóch grup da¬ 
nych 

4.1.1. Przypadek jednej zmiennej 

Przypuśćmy, że mamy do czynienia z dwoma populacjami nazywanymi również kla¬ 
sami, oznaczanymi umownie C 1 i C 2. Każdy osobnik tej populacji jest charaktery¬ 
zowany tylko jedną cechą nazywaną X. Rozkład tej cechy w obu populacjach jest 
inny, i generalnie jest opisany pewną funkcją f(x\Cl) w populacji pierwszej i funkcją 
f(x\C2) w populacji drugiej. W przypadku, gdy cecha X jest cechą numeryczną ty¬ 
pu ciągłego, funkcja / nosi nazwę funkcji gęstości prawdopodobieństwa. Przykładowe 
rozkłady funkcji / są pokazane na rysunku 4.1. 

Funkcje gęstości prawdopodobieństwa pokazane na tym rysunku bardzo przypo¬ 
minają rozkłady normalne, charakteryzujące się tą samą wariancją er 2 ale różnymi 
wartościami oczekiwanymi fii i p, 2 - 

Przypuśćmy, że naszym celem jest diagnozowanie przynależności do grupy 1, ozna¬ 
czającej - np. dla ustalenia uwagi — osoby ze stwierdzoną chorobą niedokrwienną 
serca (ChNS). Natomiast populacja II (Klasa C2) niech oznacza osoby, u których tej 
choroby nie stwierdzono. 
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Populacja II 

(o 

Populacja I 

X W,) \ 

yu 7 

X. ' Mi 

E\ - Procent osobników / 

tu - Procent osobników 

populacji I błędnie / 

\ populacji II błędnie 

zakwalifikowanych do 

zakwalifikowanych do 

populacji 11 

populacji I 


Rysunek 4.1: Rozkłady cechy X w dwóch populacjach i błędy klasyfikacji, gdy granicą 
decyzyjną jest wartość Bp. Na osi y wykreślono gęstości p-stw. 


Niech cecha X oznacza np. poziom cholesterolu stwierdzony w surowicy krwi 
u danej osoby. Naszym celem jest diagnoza ChNS. 

Generalnie, rozkłady badanej cechy w obu badanych populacjach mogą częściowo 
nachodzić na siebie — taki przypadek pokazano na rysunku 4.1. 

Naturalną zasadą klasyfikacji jest przydzielenie (zdiagnozowanie) osobnika x do 
tej klasy k, dla której wartość /(x|C7), k = 1,2, czyli wartość funkcji gęstości wy¬ 
znaczonej dla tego osobnika, jest największa. Jest to zasada, która w statystyce nosi 
nazwę zasady największej wiarygodności. 

Jeśli obie funkcje gęstości nachodzą częściowo na siebie (populacje nie są liniowo 
separabilne), to z natury rzeczy przy diagnozie możemy popełnić pewne błędy: błąd 
Ei, że osobnika należącego do populacji 1 (klasy Cl) zaklasyfikujemy do populacji 
2, i błąd E 2 , że osobnika z klasy C2 zaklasyfikujemy do Cl. 

DOBROĆ REGUŁY KLASYFIKACYJNEJ mierzymy oczekiwaną frakcją błęd¬ 
nie zaklasyfikowanych osobników. Im ta frakcja jest mniejsza, tym reguła jest lepsza. 

4.1.2. Klasyfikacja na podstawie jednej cechy o rozkładzie nor¬ 
malnym z tą samą wariancją w obu populacjach 

Rozpatrzmy bardziej szczegółowo przypadek, gdy rozkłady w obu populacjach są 
normalne, z wartościami oczekiwanymi ji\ i p 2 i tą samą wariancją a 2 . Jest to 
przypadek przedstawiony na rysunku 4.1. 

Wyznaczmy najpierw granicę decyzyjną Bp, rozdzielającą obydwa rozkłady. Wiel¬ 
kości błędów Ei i E 2 są wtedy obliczane jako odpowiednie całki (pola pod krzywymi 
odpowiednich funkcji gęstości) począwszy od wartości Bp w lewo lub w prawo. 

Punkt x = Bp wyznaczamy jako wartość, dla której gęstości rozkładów f(x\Cl) 
i f(x\C2) są takie same. Mamy więc równość 
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1 


_ exr )\— (x ~ Ml)2 j = _1_ exv\ — ^ ~ \ 

Pi 2a 2 1 { 2a 2 1 


Logarytmując obustronnie otrzymujemy 


( x - fil ) 2 
2 cr 2 


{x- fi 2 ) 2 

2a 2 


(4.1) 


(4.2) 


Po uporządkowaniu wyrażeń, wartością x dla której ta równość jest prawdziwa, oka¬ 
zuje się być 


Mj ~ lĄ = Mi +M2 

2(fti — fi 2 ) 2 


(4.3) 


Tak więc granicą decyzyjną Bq jest tutaj średnia arytmetyczna wartości oczekiwa¬ 
nych obserwowanej wartości x w obydwu populacjach. 

Nasza decyzja, czy przy danym x zaliczyć osobnika do klasy 1 czy też drugiej 
opiera się na następującej zasadzie: 


Mi > fi 2 : Jeśli x > (fii + /J. 2 )/2, to zaliczamy osobnika do klasy 1, w przeciwnym 
przypadku do klasy 2. 

fil < fi 2 = Jeśli x < (fii + /U 2 )/2, to zaliczamy osobnika do klasy 1, w przeciwnym 
przypadku do klasy 2. 

W przypadku, gdy x = (fii + yu, 2 )/2, nasza decyzja może być oparta na jakimś zran- 
domizowanym algorytmie, czego tu nie omawiamy. 

Podkreślmy jeszcze raz, zasada ta jest słuszna tylko wtedy, gdy obydwie populacje 
mają tę samą wariancję. 


4.1.3. Klasyfikacja w przypadku rozkładów wielowymiarowych 

Przenoszą się tu zasady postępowania wprowadzone dla danych jednowymiarowych. 
Jeżeli potrafimy określić rozkład /(x|Cfc) opisujący rozkład wektora badanych cech 
(x) w populacji Ck, k = 1,2, to zaliczamy osobnika x do tej populacji, dla której 
/(x|Cfc) jest największe. 

W przypadku 2 grup osobnika x zaliczamy do populacji 1, gdy 

/(x|C' 1 ) > /(x|C 2 ). (4.4) 


Pokażemy teraz, że w przypadku gdy /(x|Cfc) przedstawia wielowymiarowy roz¬ 
kład normalny, zasadę klasyfikacji wyrażoną nierównością 4.4 możemy przedstawić 
jako warunek liniowy względem rozważanych zmiennych objaśniających. 
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Przypomnijmy najpierw, że funkcja gęstości d-wymiarowego rozkładu normalnego 
ma postać: 

/(x;r,£) = — 1^1 x exp{-i[(x-/x) T S- 1 (x-/x)]}, 

(27t) 2 “IX, j 2 z 

gdzie 

x = (xi, X 2 , ■ ■ ■, Xd) T jest wektorem danych, 

//. = (/xi, ^ 2 , ■ ■ ■, Hd) T jest wektorem wartości oczekiwanych, 

jest macierzą kowariancji między zmiennymi X\, X 2 ,..., Xd (macierz sy¬ 
metryczna, na przekątnej znajdują się wariancje tych zmiennych, a poza przekątną 
kowariancje). 

Załóżmy, że rozkład ten dla populacji 1 i populacji 2 jest charakteryzowany war¬ 
tościami oczekiwanymi fi 1 i fi 2 odpowiednio. Obydwa rozkłady mają tę samą macierz 
kowariancji. 

Naszą zasadą jest, aby osobnika x zaliczyć do populacji 1 wtedy, gdy 

/(x;/i 1 ,S) > /(x; /i 2 ,S). (4.5) 

Podstawiając do tego wzoru odpowiednie wyrażenia określające rozkład normalny 
w obu populacjach, oraz logarytmując te wyrażenia otrzymujemy: 

-^[( x -Ri) Ts_ 1 ( x -Mi)] > - ^[( x -R 2 ) Ts_ 1 ( x -/^i)]- 

Przenosząc wszystkie składniki na lewą stronę nierówności, wykonując mnożenia 
i grupując względem wektora x otrzymujemy 

(Mi - ^2) Ts ~ 1x + {-^rT s ~Vi + > 0. (4.6) 

Warunek ten możemy zapisać prościej 

w T x + wq > 0, (4.7) 

gdzie 

w T = (Mi - M 2 ) Ts_1 , ^0 =(4-8) 

Jak widać, w przypadku dwóch rozkładów normalnych granicę decyzyjną stanowi 
hiperpłaszczyzna w T x + wq = 0 wyznaczona na podstawie parametrów charaktery¬ 
zujących obywa rozkłady. 

Współczynniki w i stałą Wq można obliczyć ze wzoru (4.8) lub też otrzymać za 
pomocą sieci neuronowej o architekturze perceptronu prostego z jednym neuronem 
i liniową funkcją aktywacji. 
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Gdybyśmy mieli do czynienia tylko z dwoma zmiennymi, czyli x = (xi,x2), wte¬ 
dy granicą decyzyjną (decision boundary) rozdzielającą obydwie populacje byłaby 
prosta: 

W\X\ + W 2 X 2 + w 0 = 0 

leżąca w płaszczyźnie (xi,x 2 )- Wektor w = [wi,W2] jest wektorem ortogonalnym do 
tej prostej, a stała wq wyznacza wielkość przesunięcia tej prostej względem początku 
układu. 

Ponadto można powiedzieć, że wektor w wskazuje kierunek populacji Cl. 

Fakty te można było zaobserwować w module demonstracyjnym nnd4db decision 
boundaries. Plansza z tego modułu jest pokazana na rysunku 4.2. 



Rysunek 4.2: Prosta dyskryminacyjna dla klasyfikacji 2 grup oznaczonych czarnymi 
(I) i białymi (II) kółkami. Równanie prostej rozgraniczającej: 1.41x1 + 1.41x2 — 2.12 = 
0. Wektor w = [1.41,1.41] jest wektorem ortogonalnym do prostej rozgraniczającej 
i jest skierowany na tę stronę pólpłaszczyzny , po której znajdują się punkty populacji 
wyróżnionej, u nas: czarne punkty 


4.1.4. Prawo Bayesa i obliczanie prawdopodobieństw a poste¬ 
riori 

Poprzednio chodziło nam ogólnie o zaklasyfikowanie osobnika x do jednej z dwóch 
grup danych. Obecnie chodzi nam o coś więcej: mianowicie o określenie prawdopo¬ 
dobieństwa, że osobnik scharakteryzowany wektorem x należy do określonej (&-tej) 
grupy danych. Prawdopodobieństwo to będziemy nazywać prawdopodobieństwem a 
posteriori i oznaczać p(Ck | x). 
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Prawdopodobieństwo to jest wyznaczane według prawa Bayesa 1 . 

Prawo to opiera się na dwóch następujących zasadach obowiązujących przy okre¬ 
ślaniu tzw. prawdopodobieństwa łącznego pary zmiennych (x, C), gdzie C oznacza 
zmienną określającą podział na kategorie (u nas C przyjmuje wartości C = {Ci, C 2 }), 
a x oznacza oznacza przestrzeń obserwowanych zmiennych na podstawie których 
chcemy dokonać klasyfikacji (zmienne te są czasem nazywane zmiennymi objaśniają¬ 
cymi) . 

Zakładamy, że zmienna C, tj. klasyfikacja na klasy {Ci, C 2 } pokrywa całą prze¬ 
strzeń x. 

Rysunek 4.3 ilustruje obliczanie prawdopodobieństw całkowitych i warunkowych. 



x N 




f(x\C\) - rozkład w górnej warstwie 
f(x\C 2 ) ~ rozkład w dolnej warstwie 

f(x) - rozkład całkowity, bez podziału na warstwy klasowe 

Rysunek 4.3: Rozkład łączny pary zmiennych (x, C) i wynikające z niego rozkłady 
warunkowe 


W rachunku prawdopodobieństwa są znane dwa podstawowe prawa dotyczące 
zależności między prawdopodobieństwami całkowitymi, łącznymi i warunkowymi: 

Prawo iloczynu (Product rule): P(Cfe,x) = P(Ck | x) • p(x) 

Prawo sumowania (Sum rule): -P(Cfc, x) = p(x). 

Na podstawie tych praw zostało sformułowane tzw. prawo Bayesa: 


P(C k |x) 


P(x | C k )P(Ck ) 
P(x) 


(4.9) 


przy czym p(x) może być obliczane jako 

p(x) = J2 p (x\Ck)P(C k ). 

k 

Rev. Thomas Bayes, 1702-1761 


i 
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P-stwa występujące we wzorze Bayesa noszą następujące nazwy: 

P(Ck | x) - prawdopodobieństwo a posteriori ; mówi ono, jak wysokie jest p-stwo, 
że osobnik scharakteryzowany wektorem x został wylosowany z populacji fc-tej 
(the p-blity of an event conditional on the observations); 

p{x | Cfc) - funkcja gęstości p-stwa w populacji Ck , czyli rozkład x uwarunkowany 
kategorią Ck (conditionned on class Cl); 

P(C fc ) prawdopodobieństwo a priori , że dany osobnik pochodzi z fc-tej populacji; 
p-stwo to jest wyznaczane na podstawie innych informacji, bez uwzględniania 
wektora obserwacji x dla tego osobnika (probability specified before seeing the 
data, and so based on prior experience or belief); 

p(x) - p-stwo całkowite zmiennej x, bezwarunkowe, tzn. bez uwzględniania podziału 
na klasy; jest to p-stwo brzegowe dla rozkładu łącznego (x, C ). 

Bishop ([ ], str. 82-83, również str. 232-233) pokazał, że jeżeli rozkłady warun¬ 
kowe p(x | C k ) należą do tzw. rodziny wykładniczej, to prawidłowo wytrenowana sieć 
neuronowa podaje na wyjściu właśnie prawdopodobieństwa a posteriori P(Ck | x"): 

y" = y( w ; x") = (y",..., y"), 


gdzie 


Vk ~ P(Ck I x "), k = l,...,c. 


4.1.5. P-stwa a posteriori dla 2 klas i funkcja logistyczna 

Dla przypadku gdy liczba klas c = 2, a nas interesuje klasa C\ (np. ChNS), to — 
używając reguły Bayesa — p-stwo a posteriori należenia osobnika x do klasy Ci 
możemy zapisać następująco: 


P(C : 


i |x) = 


p(x | C\)P{C\) 


p(x | Ci)P(Ci) + p(x | C 2 )P(C 2 ) 1 + exp( - a) ’ 


(4.10) 


gdzie 


p(x|C7 1 )P((7 1 ) 
p( x | C 2 )P(C 2 )' 


(4.11) 


Łatwo pokazać (por. [1], str. 233-234), że jeżeli rozkład p(x | G{) da się przedstawić 
jako rozkład należący do rodziny wykładniczej, czyli w postaci 


p(x | Ci) = exp {A(O k ) + B(x, </>) -I- 0%x} 


to obliczoną wartość a można przekształcić do następującej 

a = w T x + w o, 
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gdzie 

W = 01-02, = A(6\) — 4.(0 2 ) + ln 

W przypadku rozkładów normalnych parametry 0i i 0 2 odpowiadają wartościom 
oczekiwanym ji 1 '£~ 1 i /x 2 S _1 odpowiednio. 

Tak więc: 

jeśli obliczone dane mają rozkłady p(x | Ci) i p(x | C 2 ) należące do rodziny rozkładów 
wykładniczych, 

to obliczenia za pomocą sieci neuronowych z 1 neuronem i logistyczną funkcji akty¬ 
wacji są równoważne obliczeniom Bayesowskiego p-stwa a posteriori. 

Obliczone p-stwo a posteriori p(C\ | x) zależy od kombinacji liniowej wartości 
zmiennych zapamiętanych w wektorze x (a = w T x) za pośrednictwem logistycznej 
funkcji aktywacji. Statystycy mówią, że jest to uogólniony model liniowy z funkcją 
pośredniczącą f(a) (link function). W naszym przypadku funkcją pośredniczącą jest 
funkcja logistyczna. 

Z powyższych rozważań wynika, iż w wielu przypadkach należy spodziewać się, 
że na wyjściu sieci będą pojawiać się wartości y n które będą — za pośrednictwem 
funkcji logistycznej - zależały od zmiennej a będącej kombinacją liniową zmien¬ 
nych wejściowych x\,...,Xd- Mówimy w takim przypadku, że mamy do czynienia 
z uogólnionym modelem liniowym (ang. generalized linear model, GLM). 


4.2. Funkcja błędu E i funkcje aktywacji 

4.2.1. Funkcja błędu wyprowadzona z zasady największej wia¬ 
rygodności 

Oznaczmy wiarygodność 

L = ]]p( t"|x") 

n 

Symbol p(-) oznacza tu prawdopodobieństwo dla zmiennych dyskretnych lub funkcję 
gęstości p-stwa dla zmiennych ciągłych. 

Generalnie p(-) zależy od jakichś parametrów. Parametry te należy tak dobrać, 
aby zmaksymizować L. 

Zamiast maksymizować L możemy minimizować logarytm wiarygodności L wzięty 
ze znakiem przeciwnym: 


E= - ln L= -^lnp(t”|x”). 


n 


(4.12) 
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Bardziej konkretne określenie funkcji E zależy od postaci funkcji p( t” | x”) ozna¬ 
czającej prawdopodobieństwo warunkowe wystąpienia zdarzenia losowego t pod wa¬ 
runkiem zdarzenia x. 

4.2.2. Przypadek klasyfikacji do 2 grup 

W tym przypadku zmienna t jest jednowymiarowa i może przyjmować tylko wartości 
1 lub 0. Wobec tego p( 11 x) możemy zapisać w postaci 

p(t I x ) = y ł (l — y) 1_t , 

gdzie y oznacza ogólnie prawdopodobieństwo wystąpienia „1”. 

Powyższy zapis obejmuje zarówno przypadek wystąpienia „1” jak i „0”. 

Dla n-tego osobnika możemy powyższe zapisać jako 

p(t n | x") = (y n )*"(l — t/") 1_ * n . 

Wiarygodność zaobserwowania ciągu {x n ,£ n }, n = 1, N wynosi wtedy 

N 

E = -'£{t n lny n + (l-t n )(l-y n )} (4.13) 

n= 1 

Kryterium to nosi nazwę entropii krzyżowej (cross-entropy). 

Generalnie wartość y n jest funkcją zmiennych objaśniających x” i dodatkowych 
parametrów opisujących przyjęty model. W przypadku obliczeń przez sztuczne sieci 
neuronowe parametrami modelu są wagi w: 

j"=j(x",w). 

W przypadku sieci jednowarstwowej z jednym neuronem i logistycznej funkcji ak¬ 
tywacji wartość y n jest funkcją zmiennej a n będącej liniową kombinacją zmiennych 
wejściowych x n . 

dE _ y n — t n 
dy n y n {\ — y n ) 

The absolute minimum is at y n = t n for all n. 

V = g(a), gdzie g jest f. logistyczną, to g'(a) = g(a){ 1 - g{a)). 

5 n dE dE dy n 

da n dy n da n 

Tak więc gradient tutaj określonej funkcji błędu wyznaczony dla wartości w = 
w( fc ) wynosi 

V£| w(fc) =5>(x";w«)-i"]x". 

n 

W pakiecie NETLAB do minimizacji funkcji błędu określonej wzorem (4.13) uży¬ 
wa się metody IRLS ( Iterative Reweighted Least Sguares). 
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4.2.3. Przypadek klasyfikacji do c grup 

W tym przypadku t = [t\,... , t c ] jest wektorem o c składowych, spełniającym re¬ 
strykcje: tk € {0,1}; "}Zk^k = 1- Funkcja wiarygodności otrzymania wektora t pod 
warunkiem zaobserwowania x wyraża się wzorem 

L=p {t | x) = X\ yl k , 

k 

gdzie = 1- Wzór ten wynika z określenia rozkładu multinomialnego. War¬ 

tości yk wyrażają prawdopodobieństwa a posteriori, że osobnik scharakteryzowany 
wektorem x ra będzie należał do fc-tej klasy. 

Ujmując to jeszcze inaczej, możemy powiedzieć że symbol yk oznacza tu ogól¬ 
nie p-stwo pojawienia się „1” w fc-tej klasie (fc = 1,... ,c). Prawdopodobieństwo to 
może być różne dla różnych osobników i zależeć w szczególności od jego wartości x, 
oraz od dodatkowych parametrów, którymi w przypadku obliczeń za pomocą sieci 
neuronowych są wagi w: 

y n = y( x ”; w) 

Jeżeli rozpatrujemy próbkę N osobników charakteryzującymi się wartościami x” 
oraz t", n = 1,..., IV, to wiarygodność otrzymania takiej próbki możemy zapisać: 

L=n^ n \* n )=nmr k . 

n n k 

Wagi w i bias bo mają być wybrane tak, aby zmaksymizować wiarygodność. 
Logarytmując i biorąc wynik logarytmowania ze znakiem ujemnym, otrzymujemy 


E =-InL =■ ln( 2 /fe)- ( 4 -15) 

n k 


4.2.4. Minimalizacja błędu 

Określone wyżej funkcje błędu dla klasyfikacji do 2 lub c > 2 klas określone wzorami 
(4.13) i (4.15) zależą od wag w które należy wyznaczyć iteracyjnie w procesie uczenia 
(treningu) sztucznej sieci neuronowej. 

Funkcja błędu E jest — poza nielicznymi przypadkami - funkcją nieliniową wek¬ 
tora wag — wobec czego nie można wag wyznaczyć bezpośrednio. Funkcja ta jest 
czasami nazywana funkcją celu. 

Jeżeli E jest różniczkowalna, to możemy stosować jakieś metody gradientowe. 
Zasada jest taka, że dysponując jakimś wektorem w(fc) otrzymanym w fc-tyrn kroku, 
„poprawiamy” ten wektor poruszając się o małą odległość w kierunku, w którym 
funkcja E maleje najbardziej, czyli w kierunku ujemnego gradientu, tj. w kierunku 
W E. 



4.3. Procedury klasyfikacyjne w pakiecie Netlab 
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W ten sposób w kolejnych krokach k mamy nadzieję zbliżać się do minimum 
funkcji E. Ciąg kolejnych przybliżeń konstruujemy wg zasady: 


W ( fe+1 ) = w<*> - r 

O w 




Współczynnik nosi nazwę współczynnika uczenia (learning ratę). 

Jednak — metoda ta może zaprowadzić nas do lokalnego minimum, i w praktyce 
jest znanych i stosowanych wiele innych algorytmów optymalizacyjnych. Przykłado¬ 
wo, w pakiecie Netlab można używać takich metod, jak metoda sprzężonych gra¬ 
dientów (conjugate gradients), metoda skalowanych sprzężonych gradientów (scaled 
conjugate gradients), metoda ąuasi-Newton. 


4.2.5. Funkcje aktywacji używane w zagadnieniach klasyfika¬ 
cyjnych 

Generalnie używa się tu 3 typów funkcji: liniowej, logistycznej i softmax. 

Funkcja logistyczna ma postać: 

1 


gdzie a = w T x oznacza wynik sumowania bodźców przez neuron. 

Funkcja softmax jest określona w kontekście kilku wyjść, czyli warstwy wyj¬ 
ściowej składającej się z kilku neuronów sumujących bodźce wytwarzane przez po¬ 
przednią warstwę. 

Niech aj, = w^x oznacza wynik sumowania wykonany przez fc-ty neuron, k = 
1,..., c. Wtedy aktywacja k -tego neuronu za pomocą funkcji softmax jest określona 
następująco: 


Vk 


Efc e afc ' 


(4.16) 


Jak widzimy, funkcja softmaz zapewnia nam, że 2/fc — 1> co oznacza, że wyniki 
■yk mogą być interpretowane jako prawdopodobieństwa a posteriori. 


4.3. Procedury klasyfikacyjne w pakiecie Netlab 

Pakiet Netlab — zaimplementowany jako zestaw funkcji w MatLabie w postaci tzw. 
m-files posiada implementację kilku modeli klasyfikacyjnych: GLM, MLP, KNN 
i GMM. Modele te (za wyjątkiem KNN) są przetwarzane jako pewne specyficzne 
struktury danych. Ułatwia to dostęp do odpowiednich informacji i zapobiega pod¬ 
stawianiu przez użytkownika modelu niewłaściwych funkcji. 

Każdy z wymienionych modeli (za wyjątkiem KNN) posiada funkcję typu con- 
structor, który tworzy i inicjuje odpowiednią dla tego modelu strukturę danych 
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(obiekt) o nazwie net z odpowiednim trzy-literowym przedrostkiem (np. glmnet, mlp- 
net, gmmnet) który używany jest we wielu towarzyszących modelowi funkcjach. 

Praca modelu, a przede wszystkim otrzymywane jako rezultat tej pracy wyniki, 
są uwarunkowane tzw. opcjami. Opcje te są zapamiętywane w tablicy options(l:18). 
W systemie Netlab jest zainicjowana tablica foptions(l:18) o wartościach domyślnych. 
Użytkownik powinien skopiować tę tablicę do swojej tablicy options, a następnie 
już na gruncie swoich 'options’ — nadać poszczególnym elementom odpowiednie 
wartości. Najczęściej wykorzystywane są następujące opcje dotyczące procesu uczenia 
(train): 

options(l) - wartość 1 oznacza drukowanie błędu po wykonaniu każdej iteracji; 

options(2) - wymagana dokładność wag; 

options(3) - wymagana dokładność funkcji błędu E ; 

options(14) - maksymalna liczba iteracji. 

4.3.1. Metoda GLM — czyli uogólnionego modelu liniowego 

Model ten opisuje sieć jednowarstwową; warstwa ta jest jednocześnie warstwą wyj¬ 
ściową (warstwa wejściowa nie liczy się). Tworzenie struktury: 

net = glm(n_in, n_out, actfn) lub 

net = glm(input_dim, n_out, actfn, prior) 

Struktura 'net' przy modelu GLM: 


type 

'glin’ 


n_in 

liczba wejść (liczba cech danych) 


n_out 

liczba wyjść sieci 


nwts 

liczba wszystkich współczynników wagowych i biasów 

actfn 

f. aktywacji: string ’linear’, ’sigmoid’, 

’softmax’ 

wl, bl 

tablica wag i biasów występujących w 

pierwszej warstwie (jedynej) 


W przypadku używania złożonego modelu Bayesowskiego struktura glmnet prze¬ 
widuje dalsze pola określające parametry Bayesowskie: pola prior i beta. 

Dalsze czynności to trenowanie sieci i sprawdzenie jej skuteczności (tj. możliwości 
klasyfikacji tych samych lub innych danych testowych): 

[net, options]= glmtrain(net, options, data, targets) 

[y, a] = glmfwd(net, testdata); 

Przykładowy skrypt wykonujący klasyfikację za pomocą metody GLM: 
net=glm(3,1, 'logisticO ; 

data=w65(: , [2 3 4]); °/ 0 bierzemy tylko kolumny 2, 3, 4 

target=w65(:,5); °/ 0 bierzemy kolumnę 5 tablicy z~danymi 

options=zeros(l,18); options(1)=1; options(14)=30; 
net=glmtrain(net,options,data, target); 
tdata=gr7 ( : , [3 5 6] ) ; °/ 0 dane testowe 

y=glmfwd(net, tdata); °/ 0 wyniki sieci dla danych testowych 



4.3. Procedury klasyfikacyjne w pakiecie Netlab 


47 


4.3.2. Metoda MLP — dla dwuwarstwowego perceptronu 

Implementacja w Netlabie uwzględnia tylko perceptrony dwuwarstwowe, tj. zawiera¬ 
jące jedną warstwę ukrytą i jedną warstwę wyjściową. Tworzenie struktury: 

net = mlp(n_in, n_hidden, n_out, actfn) 
lub też — w przypadku używania złożonego modelu Baysowskiego: 
net = mlp(input_dim, n_hidden, n_out, actfn, prior.beta) 

Struktura ’net’ przy modelu MLP: 


type 

’mlp’ 

n_in 

liczba wejść (liczba cech danych) 

n_hidden 

1. neuronów w warstwie ukrytej 

n_out 

liczba wyjść sieci 

nwts 

liczba wszystkich współczynników wagowych i biasów 

actfn 

f. aktywacji: string dinear, ’sigmoid’, ’softmax’ 

wl, bl,w2,b2 

tablica wag i biasów występujących w obu warstwach 


Struktura może zawierać dalsze pola określające parametry złożonego modelu 
Baysowskiego. 

Dalsze czynności to trenowanie sieci i sprawdzenie jej skuteczności (tj. możliwości 
klasyfikacji): 

[net, options, varargout]= netopt(net, options, x, t, alg); 

[e, edata, eprior] = mlpfwd(net, testdata); 

Przykładowy skrypt wykonujący klasyfikację za pomocą metody MLP: 

nh=5; °/* 1. neuronów w~warstwie ukrytej 

net=mlp(3,nh,1,’logistic’); 
data=w65(:,[2 34]); target=w65(:,5); 
options=zeros(l,18); options(1)=1; options(14)=30; 

[net, options] =netopt(net, options, data, target, 'conjgradO; 

7, innym optymizerem może byc np. : ’scg ł lub 'ąuasinew’ 
tdata=gr7(:, [3 5 6]); 

options5=options; 7« to jest na ogol niepotrzebne 

y_mlp5=mlpfwd(net,tdata); 

7, tab= [gr7( : , 1) y_mpll y_mlp3 y_mlp5] ; 


4.3.3. Metoda KNN — najbliższych K sąsiadów 

4.3.4. Metoda GMM — bazująca na mieszaninie rozkładów 




Rozdział 5 


Algorytmy minimalizacji 
błędu 

5.1. Znajdywanie minimum funkcji błędu 

Funkcja błędu E = E( w) jest w ogólnym przypadku nieliniową funkcją wektora wag 
w = [uą, iv 2 ,..., wg], gdzie S jest sumaryczną liczbą wag łącznie z wartościami pro¬ 
gowymi (biasami). Funkcja E przedstawia najczęściej sumę kwadratów różnic między 
wartościami docelowymi (t n ) a wartościami y”(x n ; w) wyprodukowanymi przez sieć; 
czasami — szczególnie w zagadnieniach klasyfikacyjnych — jest to ujemny logarytm 
z prawdopodobieństw aposteriori lub z wiarygodności P(t n \bfx n ) otrzymania warto¬ 
ści docelowych t n pod warunkiem zaobserwowania wektora zmiennych objaśniających 


5.2. Klasy algorytmów minimalizacyjnych 

Generalnie możemy tu wyróżnić trzy klasy: 

1. Algorytmy klasyczne. Są to algorytmy związane z zagadnieniami minimaliza- 
cyjnymi w ogóle. Nie wykorzystują specyfiki sieci neuronowych i są stosowane 
w wielu innych zagadnieniach optymalizacyjnych. Najczęściej są to algorytmy 
gradientowe, chociaż są znane również algorytmy bezgradientowe. 

2. Algorytmy związane specyficznie z sieciami neuronowymi. Zalicza się tu przede 
wszystkim algorytm tzw. propagacji wstecznej czyli backpropagation. Algorytm 
backpropagation bierze pod uwagę rozkład neuronów (a tym samym przypisa¬ 
nych im wag) w poszczególnych warstwach. 



50 


Rozdział 5. Algorytmy minimalizacji błędu 


3. Algorytmy z elementami stochastycznymi. Należą tu w pierwszym rzędzie: algo¬ 
rytm symulowanego wyżarzania, algorytmy genetyczne i algorytmy ewolucyjne. 

W dalszym ciągu omówimy jedynie kilka algorytmów klasycznych wykorzystują¬ 
cych metody gradientowe. 


5.3. Algorytmy klasyczne służące minimalizacji funk¬ 
cji E(w) 

Generalnie rzecz biorąc, są to algorytmy iteracyjne. 


5.3.1. Zasada algorytmów iteracyjnych 

Z jakiegoś początkowego przybliżenia wektora wJ°) konstruujemy w kolejnych kro¬ 
kach (iteracjach) ciąg kolejnych przybliżeń wektora w który ma dać minimum funkcji 
błędu E{ w). Przypuśćmy, że po k krokach otrzymaliśmy wartość W następnym 
kroku, (k + l)-szym, konstruujemy następne przybliżenie, (k + l)-sze, według zasady: 

w (‘+ 1 )=w<‘)+Aw( t ). (5.1) 

Cała sztuka algorytmu leży w tym, aby w skuteczny sposób znaleźć poprawkę 
Poprawka ta polega na ogól na znalezieniu pewnego wektora p( fc ) i dodaniu 
części tego wektora do posiadanego aktualnie wektora (nie możemy dodać za 
dużo, bo moglibyśmy „przeskoczyć” właściwe minimum). Tak więc podstawiamy 


Aw< fc >=łjp< fc >. (5.2) 

Dalej poznamy różne metody znajdywania kierunku p. W najprostszym przy¬ 
padku może to być ujemny gradient funkcji E w punkcie ~w^ k \ co oznacza że 
podstawiamy p( fc ) = —X7E\ w (k). Wtedy poprawka Aw^) przyjmie postać: 


Aw (fc > 


—r)VE 




(5.3) 


Współczynnik rj nosi nazwę współczynnika uczenia. Powinien on przyjmować 
wartości z przedziału [0, 1]. Na ogół wprowadza się go jako funkcję malejącą 
numeru iteracji. 

Jeżeli wartości ij są zbyt małe, to postęp w zbliżaniu się do minimum jest zbyt 
wolny; jeżeli 77 jest zbyt duże, to wartości w mogą zmieniać się zbyt chaotycznie. 
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• Aby przeciwdziałać zbyt chaotycznym zmianom, wprowadza się do wzoru na 
dodatkowy człon, tzw. momentum lub pęd: 


Aw (fc) =-r]VE +p(w (fe) - w (fc - 1} ) 



(5.4) 


momentum 


Dodatkowy człon momentum zapobiega raptownym zmianom kolejno wyzna¬ 
czonych wartości wW; można powiedzieć, że do pewnego stopnia zachowuje 
kierunek (tendencję) zmian. 

Określenie odpowiednich wartości rj i p jest sprawą wyczucia i doświadczenia. 
Technika bold driuer techniąue ([1], str. 269), pozwala wyznaczyć jednocześnie 
współdziałające wartości rj i p. 


5.3.2. Aproksymacja funkcji błędu za pomocą rozwinięcia w sze¬ 
reg Taylora — i co z tego wynika 

Rozwinięcie funkcji E(yv + p) dookoła punktu w: 


E(yt + p) = E( w) + g(w) T p + -p 7 H(w)p + 0(h 3 


(5.5) 


gdzie: 


g(w) = VE 



' dE 

dE 


W 

dwi 

w”"’ dw s 

W- 


H = 


d 2 E 

dwidwj 


> i,j = l,...,S. 

W 


(5.6) 

(5.7) 


Wyrażenie g(w) nazywa się gradientem funkcji (powierzchni) E w punkcie w i określa 
płaszczyznę styczną do tej powierzchni. 

Macierz H nazywa się hesjanem i zawiera w sobie informacje o krzywiźnie po¬ 
wierzchni. 


Różne algorytmy służące minimizacji w dużym stopniu korzystają z przedstawio¬ 
nego wyżej rozwinięcia (5.5); przy czym niektóre odmiany tych algorytmów biorą 
tylko człony liniowe, a inne zarówno człony liniowe jak i kwadratowe. 


5.3.3. Algorytm największego spadku 

W rozwinięciu (5.5) ograniczamy się tylko do członu liniowego rozwinięcia: 


E( w (fe) + p (fc) ) = E{ w (fe) ) + g(w (/c) ) T p (fc) + 0{h 2 ). 


( 5 . 8 ) 

















52 


Rozdział 5. Algorytmy minimalizacji błędu 


Chcemy, żeby 

E(w^ k+1 ^) < E( w<‘>). 

Będzie tak, gdy 

g(w«fp«<0. 

Wynika stąd, że będzie tak (tj. nastąpi obniżenie błędu), gdy przyjmiemy: 

P (fe) = -g(w (fc) ). 

W ten sposób otrzymaliśmy kierunek poprawki. Samej poprawki Aw^ dokonu¬ 
jemy korzystając z wzoru 5.3 lub 5.4. 

Wadą metody jest niewykorzystywanie informacji o krzywiźnie funkcji (powierzch¬ 
ni błędu). Ponadto, w okolicy minimum gradient bywa mały, a więc zbliżanie się do 
minimum staje się bardzo powolne. 

Zaletą metody jest jej prostota, małe wymagania co do pamięci, stosunkowo 
mała złożoność obliczeniowa. 

Ponadto — można uzyskać poprawę efektywności metody przez zastosowanie mo- 
mentum. 

Dalsze modyfikacje tej metody to: reguła delta-bar-delta (wyznacza m.in. w spe¬ 
cjalny sposób współczynnik uczenia rf) i metoda ąuickprop Fahlmana. 

5.3.4. Metoda Newtona 

Korzystamy tutaj z przybliżenia kwadratowego funkcji E( w): 

E(w (k) + p (fc) ) w E{ w (fe) ) + g(w (fe) ) T p (fe) + ^(p (fe) ) T H(w (fe) )p (fe) . (5.9) 

Szukamy takiego wektora p( fc ), który by — dodany do wektora — pozwolił nam 
osiągnąć minimum powierzchni błędu. 

Jeżeli osiągnęliśmy minimum, to pochodna w tym punkcie powinna być równa 
zero: 

dE( w (t 4p w ) 

OP “ ' 

Różniczkując obustronnie równość (5.9) ze względu na wektor p^, i korzystając z te¬ 
go, że pochodna lewej strony wyznaczana w minimum jest równa zeru, otrzymujemy: 


0 = g (w (fc) ) +H(w (fc) )p (fe) . 


Wynika stąd, że jako poszukiwany wektor kierunkowy p( fc ) doprowadzający nas do 
minimum należy przyjąć 
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Aby zapewnić ten warunek, należy w każdym cyklu określić wartość gradientu g 
oraz Hesjanu H w punkcie znanego (ostatniego) rozwiązania w^ k \ 

Wtedy: 

w (k+i) = w (fc) _ [H(wW)]- 1 g(wW). (5.11) 

Wzór ten jest trudny do zrealizowania w praktyce, gdyż 

• wymaga w każdym kroku iteracyjnym (k) obliczenia wartości hesjanu H w punk¬ 
cie 

• wymaga odwrócenia tego hesjanu, do czego jest potrzebna dodatnia określoność 
tego hesjanu i to w każdym kroku iteracji - czego sam algorytm nie jest 
nam w stanie zapewnić. 

W tym stanie rzeczy w praktycznych implementacjach rezygnuje się z wyznacza¬ 
nia dokładnej wartości hesjanu, a zamiast tego korzysta się z jego wartości przybli¬ 
żonej, która w dalszym ciągu będzie oznaczana jako H(w^ fc ^). 

Jedną z najpopularniejszych metod służących temu celowi jest metoda zmiennej 
metryki. 

5.3.5. Metoda zmiennej metryki — ąuasi-Newton 

W metodzie tej w każdym kroku iteracji modyfikuje się hesjan H z kroku poprzed¬ 
niego o pewną poprawkę. 

Poprawka może być tak dobrana, aby aktualna wartość hesjanu H(w^) przybli¬ 
żała krzywiznę funkcji celu E zgodnie z zależnością: 

H(w (fc) )(w (fc) - w (fe-1) ) = g(w (fc) ) - g(w (fe - 1) ), 
gdzie g(-) oznacza odpowiedni gradient. 

Na podstawie powyższego założenia można otrzymać wzory określające hesjan 
w kroku /c-tym w zależności od wartości hesjanu w kroku ( k — l)-szym oraz od przy¬ 
rostów gradientu i optymizowanej zmiennej w. Przybliżenie hesjanu H(w^ fc ^) powinno 
być tak skonstruowane, aby zapewnić istnienie odwrotności V = [H(w( fc ))] _1 . 

Osowski [12] (str. 57) podaje dwa algorytmy posiadające tę własność. Są to: 

- algorytm Broydena-Goldfarba-Fletchera-Shanno (BFGS) oraz 

- algorytm Dawidona-Fletchera-Powella (DFP). 

Metoda zmiennej metryki charakteryzuje się zbieżnością superliniową; jest więc 
znacznie lepsza od liniowo zbieżnej metody największego spadku. 

Według Osowskiego (str. 57) metoda ta jest obecnie uważana za jedną z najlep¬ 
szych metod optymalizacji funkcji wielu zmiennych. 
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Jej wadą jest stosunkowo duża złożoność obliczeniowa wynikająca z konieczności 
wyznaczania S 2 elementów hesjanu, a także duże wymagania co do pamięci przy 
przechowywaniu macierzy hesjanu. Jej skuteczność sprawdzono na komputerach oso¬ 
bistych dla sieci nie przekraczających tysiąca parametrów. 

5.3.6. Algorytm Levenberga-Marquardta 

Jest inną, bardzo popularną odmianą newtonowskiej metody optymalizacji. War¬ 
tość dokładna hesjanu ze wzoru (5.10) zostaje zastąpiona wartością aproksymowa- 
ną H(w( fe )) określaną na podstawie informacji zawartej w gradiencie, a także przy 
uwzględnieniu czynnika regularyzującego (Osowski, str. 58-59). 

5.3.7. Metoda sprzężonych gradientów 

Można by usiłować budować ciąg kolejnych kierunków określających kolejne po¬ 
prawki Aw( fc ) (por. wzór (5.2)) jako ciąg wektorów pO) wzajemnie ortogonalnych. 
Zapewnia to przeszukiwanie przestrzeni optymalizowanych parametrów w niezależ¬ 
nych kierunkach. Jednak w takim przypadku nie wykorzystujemy informacji o krzy- 
wiźnie powierzchni błędu. Aby to zrobić, konstruujemy ciąg kolejnych przybliżeń, 
czyli ciąg wektorów p( fc ) w ten sposób, żeby one były wzajemnie sprzężone względem 
hesjanu H lub jego przybliżenia. 

Definicja. Wektory p j i p, są wzajemnie sprzężone względem macierzy G, jeśli 
jest spełniony warunek: 

pjGp i = 0. 

Osowski (str. 59) pokazuje, że kolejny kierunek p- fe ) można otrzymać jako kom¬ 
binację liniową aktualnego gradientu g^ fe) i kierunków p (0 ), ..., pA- 1 ) wyznaczonych 
w poprzednich iteracjach: 


k-1 

P (fe) = -g k + ^2f3kjP b) , 
j =o 

gdzie g k jest gradientem funkcji błędu wyznaczonym w punkcie w 

Po uwzględnieniu warunków ortogonalności i warunków sprzężenia względem he¬ 
sjanu H otrzymujemy 

P (fc) = -gfc + )9fc- 1 p (fc - 1) , 

gdzie współczynnik f3k -1 może być wyznaczony np. według jednej z następujących 
reguł (Osowski, str. 60, Bishop, str. 280-281): 


Pk-l 


gfc (gfc ~ gfc-i) 
sl-iSk-i 


(Polak-Ribiere) 
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R _ gfc (gfc ~ gfc-i) 

Pk-i = , — Sk ~^ , (Hestenes-Stiefel) 

(P fc_1 ) T (gfc - gfc-i) 

T 

(3k-i = r k ^ k — (Fletcher-Reeves) 

§fe—igfc—i 


Tak więc — korzystając z zasady sprzężenia — możemy konstruować kie¬ 
runki poszukiwań. 

Następnym etapem jest znalezienie minimum funkcji na danym kierunku p A 
przez tzw. minimalizację kierunkową. Wykorzystujemy tu aproksymację kwadratową 
wynikającą z rozwinięcia Taylora (5.9) postaci: 


%W+a P W ) w £(w( fe ))+g k T (wW +np (fc) ) + ^(w<*> +ap (fc) ) T H(wW +ap (fc >). 

Teraz szukamy takiej wartości a która dałaby minimum funkcji błędu na kierunku 

( w ( fe ) + ap( fc )). 

Obliczając pochodną cząstkową względem a i przyrównując ją do zera otrzymamy 

(fc) = P (fc) gfc (512) 

min ( p (fc) )T H p(0 ' 1 ’ 

Jak można zauważyć, do wyznaczenia wielkości a^ln wystarczy znajomość HpA') — 
który można otrzymać pośrednio, bez oddzielnego wyznaczania całego hesjanu H. 

Typowo, obydwa te kroki (wyznaczanie nowego, sprzężonego kierunku i szu¬ 
kanie stałej o:fc dającej minimum na tym kierunku) są stosowane naprzemiennie. 

Własności metody. Według Osowskiego (str. 60) metoda sprzężonych gradientów 
wykazuje zbieżność zbliżoną do liniowej i z tego powodu jest mniej skuteczna niż 
metoda zmiennej metryki, ale zdecydowanie szybsza niż metoda największego spad¬ 
ku. Stosuje się ją powszechnie jako jedyny skuteczny sposób algorytmu optymalizacji 
przy bardzo dużej liczbie zmiennych sięgających nawet kilkudziesięciu tysięcy. 

5.3.8. Metoda gradientów sprzężonych z regularyzacją (scaled 
conjugate gradients) 

Algorytm scg jest dość skomplikowany. Jest on opisany u Osowskiego na str. 68-71 
i Bishopa str. 282-285, jednak opisy w obu tych źródłach różnią się. 

Generalnie, algorytm ten wiąże się z nazwiskiem Moliera i pozwala unikać bez¬ 
pośredniego szukania minimum na danym kierunku przez tzw. linę search. Algorytm 
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stosuje regularyzację dla uniknięcia trudności wynikających z faktu, że hesjan rozwa¬ 
żany w tej metodzie (wynikający z aproksymacji kwadratowej funkcji błędu w punkcie 
w( fe i) na ogół bywa niedodatnio określony lub źle uwarunkowany. 

Regularyzacja jakiejś wielkości wiąże się generalnie z dodaniem pewnego warunku 
na tę wielkość. Z dodaniem nowego warunku wiąże się dodanie nowego parametru do 
rozpatrywanego modelu. 

W omawianym zagadnieniu regularyzacja może zawierać dwa momenty: 


i) uczynić hesjan dodatnio określony przez dodanie do macierzy H pewnej krotno¬ 
ści A macierzy jednostkowej; tym samym w dalszym ciągu będzie rozpatrywana 
macierz H +Al; 


ii) 


wprowadzić regularyzację do formuły na długość kroku a^ wyznaczanego we¬ 
dług wzoru (5.12), co powoduje, że krok ten będzie wyznaczany z następującego 
warunku 


(k) = _ P (fc) gk _ 

(p( fc )) T Hp( fe ) +Afc||pi fc i 


(5.13) 


Parametr A& występujący w warunku regularyzacji jest dobierany indywidualnie 
w każdym kroku iteracji. 


Jak pisze Bishop (str. 285), rezultaty symulacji komputerowych wykazują, że w nie¬ 
których przypadkach algorytm scg może dać znaczącą poprawę w szybkości obliczeń 
w stosunku do zwykłego algorytmu sprzężonych gradientów. 






Rozdział 6 


Sieci o radialnych funkcjach 
bazowych czyli sieci typu 


RBF 


Sieci te są opisane np. w książkach Osowskiego ([12] rozdział 5) i Bishopa ([ ] również 
rozdział 5). 

6.1. Architektura i funkcjonowanie sieci typu RBF 

Sieci tego typu służą najczęściej do aproksymacji zmiennej numerycznej. Są używane 
również w zagadnieniach klasyfikacyjnych. Z różnych twierdzeń (np. Poggio i Giro- 
siego) wynika, że taka sieć może służyć jako uniwersalny aproksymator i przybliżyć 
dowolną funkcję ciągła z dowolną dokładnością. 

Sieci RBF są sieciami dwuwarstwowymi, z przepływem informacji do przodu, bez 
wymiany informacji między neuronami w jednej warstwie. 

Tak więc sieć typu RBF oprócz warstwy wejściowej dostarczającej dane x składa 
się: 

1. z warstwy H neuronów ukrytych, 

2. z warstwy wyjściowej liczącej ogólnie K wyjść. 

W poprzednio omawianych sieciach jedno- i wielowarstwowych neurony występu¬ 
jące w liczących się warstwach spełniały rolę sumatorów impulsów (bodźców) docho¬ 
dzących do nich z poprzedniej warstwy. Rezultat sumowania byl następnie przetwa¬ 
rzany przez tzw. funkcję aktywacji i powodował pobudzenie neuronu na określonym 
poziomie, co z kolei miało ten skutek, że dany neuron przekazywał swoją aktywację 
jako bodziec numeryczny — do następnej warstwy, lub na wyjście. 
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Rozdział 6. Sieci o radialnych funkcjach bazowych (typu RBF) 


W przypadku sieci typu RBF każdy neuron warstwy ukrytej odgrywa rolę małego 
kalkulatorka: wyznaczanie aktywacji odbywa się tu na zasadzie wyznaczenia warto¬ 
ści tzw. radialnej funkcji bazowej — o której zaraz powiemy coś więcej. Poza tym, 
wszystko odbywa się tak samo, jak w zwykłych sieciach dwuwarstwowych. 

Definicja. Radialną funkcją bazową (typu RBF) nazywamy funkcję G(-) postaci 

G(x;c) = G(r(x,c)), gdzie r(x,c) = ||x - = {(x - c) T (x - c)} 1/2 . (6.1) 

Oznacza to, że wartości funkcji dla danego argumentu x - zależą tylko od odle¬ 
głości jej argumentu x od centrum c będącym parametrem tej funkcji. Czasami taka 
pojedyncza funkcja radialna jest nazywana jądrem (kernel), a parametr a szerokością 
(width) jądra. 

Stosunkowo często używaną funkcją radialną jest funkcja Gaussa 

G(x;c,cr 2 ) = exp{J^-^-}. (6.2) 

Jest to funkcja o kształcie dzwonu. Szerokość dzwonu reguluje parametr a 2 nazywany 
również parametrem gładkości. 

Przykłady innych funkcji radialnych (wraz z wykresami) można znaleźć np. w książ¬ 
ce Osowskiego, str. 168, Rys. 5.4. 

Tak więc neurony pierwszej warstwy obliczają — na podstawie podanego na wej¬ 
ściu wektora x — swoje aktywacje jako wartości 

Gi(x) = G(x,ci), G 2 (x) = G(x,c 2 ), ... , G ff (x) = G(x,cjj). 

Subskrypt u dołu symbolu G oznacza tutaj numer centrum. Tak więc symbol Gh 
oznacza funkcję radialną obliczaną względem centrum c^. 

Obliczone w ten sposób wartości Gi, G 2 ,..., Gh służą jako dane wejściowe dla 
warstwy drugiej, wyjściowej, która oblicza z nich ważoną sumę [k = 1,..., K, gdzie 
K oznacza liczbę wyjść): 

J/fc(x) = w k o + tUfciGi(x) + tUfc 2 G 2 (x) H-hWfcj?Gij(x). (6.3) 

co można zapisać w postaci skróconej: 

Uk = Wko + Wk\G\ + rUfc 2 G 2 + • • • + WkHGn- (6-4) 

Przy danych centrach ci,...,c h i parametrze gładkości er 2 pozostają nam do 
wyznaczenia tylko wektory wag wi,..., wjf używane przez warstwę wyjściową. 

Wektory te otrzymuje się jako rezultat uczenia sieci za pomocą tzw. próbki uczą¬ 
cej. Próbka taka zawiera zbiór wzorców danych przez pary |x”,t”}, dla których 
przy każdym wektorze x" znana jest prawidłowa odpowiedź t", jakiej sieć powinna 
dostarczyć. 

W dalszym ciągu omówimy bardziej szczegółowo tylko przypadek sieci typu RBF 
z jednym wyjściem. 



6.2. Wyznaczanie centrów i parametrów gładkości 
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6.2. Wyznaczanie centrów i parametrów gładkości 

Centra — w zadeklarowanej liczbie H — są rozstawiane w przestrzeni zmiennych 
objaśniających R d . Wyznaczanie centrów może się odbywać w zasadzie na trzy spo¬ 
soby: 

1. w sposób losowy; 

2. przez podział danych punktów (określonych zmiennymi objaśniającymi) na H 
klasterów; 

3. zestartowanie z liczbą klasterów równą liczebności próbki uczącej — każdy 
punkt tej próbki stanowi wtedy centrum; a następnie przez stopniową redukcję 
liczby centrów aż zostanie ich tylko H. 

Parametr er 2 wyznacza się na ogół metodą cross-validation (fc-fold cross-validation). 
Możliwy jest również algorytm który przyjmuje, że punkty próbki uczącej są mie¬ 
szaniną rozkładów gaussowskich o indywidualnych centrach i indywidualnych szero¬ 
kościach jąder. Jednoczesne wyznaczanie centrów c^, (h = 1,..., H) i odpowiadają¬ 
cych im szerokości ah jest możliwe dzięki naprzemiennemu zastosowaniu algorytmu 
EM. Algorytm takiego wyznaczania klasterów jest zrealizowany w pakiecie Netlab. 

6.3. Sieć typu RBF z jednym wyjściem i jej uczenie 

Zakładamy, że w warstwie ukrytej znajduje się H neuronów, z których każdy (np. 
neuron h- ty) ma przyporządkowane swoje indywidualne centrum określające strefę 
wpływów tego neuronu. 

Zakładamy również, że aktywacja każdego neuronu będzie się odbywać za pośred¬ 
nictwem radialnej funkcji Gaussa o parametrze a 2 wspólnym dla wszystkich neuro¬ 
nów. 

Tak więc mamy następującą sytuację 

Neuron Funkcja aktywacji 

Neuron 1 : G(||x — C!|| 2 , er 2 ) 

Neuron 2 : G(||x — c 2 1 | 2 , er 2 ) 

Neuron H : G(||x — Cflj| 2 , a 2 ) 

Możliwe są również modele, w których każdy neuron posiada swój własny para¬ 
metr a 2 . 

6.3.1. Uczenie sieci 

Ponieważ mamy tylko jedno wyjście, więc zbiór próbek uczących (wzorców) będzie 
miał postać 

{x™, t n }, n=l,...,N. 
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Rozdział 6. Sieci o radialnych funkcjach bazowych (typu RBF) 


Tym samym, dla n-tego wzorca warstwa ukryta dostarcza wartości 
Gi(x"), G 2 (x"), ... ,G ff (x"), 
z których jest obliczany wynik (odpowiedź) sieci na sygnał x": 

y(x”) = wiGi + w 2 G 2 + ■■■ w H G h + w 0 . 

Porównując odpowiedzi sieci y(x ra ) z wartościami docelowymi t n określamy ogólny 
błąd odpowiedzi E. Błąd ten jest wyznaczany najczęściej na zasadzie najmniejszych 
kwadratów: 


N N H 

E = 2 - 2/( x ")] 2 = 2 E^" - E^ G i( x ")] 2 - (6-5) 

n =1 n—1 j—0 

gdzie symbolem Go = Go(x”) oznaczono funkcję przyjmującą tożsamościowo wartość 
równą 1. 

Ponieważ wartości Gj(x") są znane, błąd E jest kwadratową funkcją wag w = 
[u>o, W\,W 2 , , Wh]- W takim przypadku (tj. gdy pod kwadratem występuje liniowa 

funkcja szukanych parametrów wo,Wi, ... ,wh), wektor wag w dający minimum 
błędu może być znaleziony explicite jako rozwiązanie układu równań liniowych (np. 
tzw. układu równań normalnych). Nie potrzeba już metod iteracyjnych — co 
czyni sieć typu RBF bardzo atrakcyjną z punktu widzenia obliczeniowego. 

Jednak - - przy rozwiązywaniu układu równań liniowych można również napotkać 
na trudności, np. otrzymany układ może być źle uwarunkowany. 

W książce Osowskiego można znaleźć — w rozdziale 5 poświęconym sieciom RBF 
kilka metod rozwiązywania układu równań liniowych wyznaczających wagi uży¬ 
wane przez takie sieci. Może to być np. pseudoodwrotność lub SYD. 


6.4. Implementacja sieci radialnych w pakiecie Ne- 
tlab 

Pakiet Netlab posiada swój własny konstruktor dla sieci RBF. Jego działanie poka¬ 
zuje moduł demonstracyjny DEMRBF1. Moduł ten pokazuje, jak można otrzymać 
przybliżenie sinusoidy za pomocą różnego typu funkcji radialnych: gaussowskiej, thin 
piąte spline (’TPS’) (G(r) = r 2 logr) oraz ’r41ogr’. 

Przed przystąpieniem do obliczeń dane powinny być wystandaryzowane. 

A oto przykładowy skrypt wykonujący obliczenia dla 6-u różnych liczebności war¬ 
stwy ukrytej. 

"/, Demonstrates calculations of a~RBF approximation of a~sine function 
‘/, Generate the matrix of inputs x and targets t. 



6.4. Implementacja sieci radialnych w pakiecie Netlab 
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randn(’State’, 42); rand(’State’, 42); 
ndata = 20; 7. Number of data points. 

noise =0.2; % Standard deviation of noise distribution. 

x = (linspace(0, 1, ndata))’; 
t = sin(2*pi*x) + noise*randn(ndata, 1); 
mu = mean(x); sigma = std(x); 

tr_in = (x - mu) .// (sigma); l Standardized data 

H= [2 4 7 12 15 18]; 7, Number of hidden neurons 
for hh=H 

7. Set up network parameters. 

nin = 1; 7. Number of inputs. 

nhidden = hh; 7. Number of hidden units. 

nout = 1; 7. Number of outputs. 

7. Create and initialize network weight and parameter vectors. 
net = rbf(nin, nhidden, nout, ’gaussian’); 

7. Use fast training method 
options = foptions; 

options(l) =0; 7. No Display of EM when making clusters 

options(14) = 10; 7. number of iterations of EM 

net = rbftrain(net, options, tr_in, t); 
y = rbffwd(net, inputvals); 

disp([’RBF training errors are for ’,... 

'nhidden H= ’, num2str(nhidden),... 

’ RBF Gaussian ’, num2str(rbferr(net, tr_in, t))]); 


end 


W rezultacie otrzymujemy następujące wartości błędu E: 


RBF training errors are for nhidden H=2 
Gaussian 0.5087 TPS 4.5182 R41ogr 4.8817 

RBF training errors are for nhidden H=4 
Gaussian 0.24256 TPS 0.70335 R41ogr 0.9769 

RBF training errors are for nhidden H=7 
Gaussian 0.14586 TPS 0.14092 R41ogr 0.26191 

RBF training errors are for nhidden H=12 
Gaussian 0.10797 TPS 0.089096 R41ogr 0.097692 
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Rozdział 6. Sieci o radialnych funkcjach bazowych (typu RBF) 


RBF training errors are for nhidden H=15 
Gaussian 0.10797 TPS 0.0032335 R41ogr 0.0073178 

RBF training errors are for nhidden H=18 
Gaussian 0.10792 TPS 0.0014216 R41ogr 0.0043882 


Wykonano również (za pomocą DEMRBF1) rysunki ilustrujące aproksymację 
funkcji sinus za pomocą sieci radialnych z H=2, 4, 7 i 15 neuronami. 




6.4. Implementacja sieci radialnych w pakiecie Netlab 
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Rysunek 6.1: Aproksymacja funkcji sinus wykonana za pomocą sieci o H=2 (góra) 
i H=4 (dół) neuronach. Sieci zostały wytrenowane na podstawie 20-elementowej za¬ 
burzonej próbki wylosowanej z sinusoidy. 
















Target Target 
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Rozdział 6. Sieci o radialnych funkcjach bazowych (typu RBF) 




Rysunek 6.2: To samo, co na rysunku 6.1, ale warstwa ukryta sieci składała się z H=7 
(góra) i H=15 (dół) neuronów. 
















Rozdział 7 


Uczenie samoorganizujące się 

7.1. Uwagi Tadeusiewicza nt. Sieci samoorganizują¬ 
cych się 

Sieci samoorganizujące się wykonują grupowanie danych. Pojawiają się dwie zasad¬ 
nicze własności powodujące istotne zmiany w funkcjonowaniu sieci: 

A. Koherencja, 

B. Kolektywność. 

Ad A. 

Zasada koherencji postuluje, żeby grupować dane w pewne klasy podobieństwa. 
Grupowanie następuje na zasadzie samoorganizacji. Proces samoorganizacji jest au¬ 
tomatyczny i spontaniczny, bez nauczyciela. Zakłada się, że komplet potrzebnych 
informacji jest zawarty w samych danych, z których jedne są podobne do siebie na¬ 
wzajem, a inne nie. 

Ad B. 

Zasada kolektywności powoduje, że to, co rozpoznaje jeden neuron, zależy w dużej 
mierze od tego, co poznają inne neurony. 

Zbiorowość neuronów (czyli ich kolektyw) może w sposób pełniejszy i bogatszy 
przetwarzać informacje, aniżeli każdy z neuronów wzięty z osobna. 

Zbiorowość odpowiednio powiązanych i współpracujących elementów stwarza moż¬ 
liwość uzyskania nowych form zachowania i nowych postaci działań znacznie bogat¬ 
szych, niż by można oczekiwać biorąc pod uwagę każdy z elementów tej zbiorowości 
z osobna. Przykłady z przyrody: Pszczoły, mrówki jako kolektyw są zdolne do 
celowych działań. 
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Rozdział 7. Uczenie samoorganizujące się 


7.2. Zasada Hebba 

7.2.1. Oznaczenia wstępne 

Rozpatrujemy prosty element sieci neuronowej 1 , w której neuron o numerze i pobiera 
w fc-tym kroku (fc-tej iteracji) dane x(fc) = x = [aą,... ,Xd] T i na tej podstawie, ko¬ 
rzystając ze swoich aktualnych wag w* = ..., Wid) T wytwarza wynik yi. Wynik 

ten jest obliczany jako 

»,=wr*=*v. 

Element sieci pokazujący jeden neuron, oznaczony numerem i, jest pokazany na ry¬ 
sunku 7.1. Zauważmy, że neuron ten posiada swój wektor wag, w. j; na podstawie 

Al 

X 2 


X 3 

X 4 



Rysunek 7.1: Model liniowego neuronu z wyjściem yi = wj x. 



w 


rzut x na w 


Rysunek 7.2: Rzutowanie wektora x na wektor w 

których oblicza swoją aktywację jako kombinację liniową XjWij, j = 1, 2,..., d, 
gdzie d jest wymiarem wektora wejściowego, czyli sygnału pre-synaptycznego; na ry¬ 
sunku wymiar ten wynosi d = 4. Otrzymane w ten sposób pobudzenie neuronu jest 
kierowane na wyjście jako tzw. sygnał post-synaptyczny yi. Tak więc, tak naprawdę, 
wartość wyjścia t/j zależy od wielkości sygnału pre-synaptycznego x oraz wektora 
wagowego wi = [uąi,... ,Wid] T przypisanego do neuronu n,. 

Zauważmy również, że gdyby wektor w., był wektorem jednostkowym, to wartość 
yi mogłaby być interpretowana jako rzut wektora x na wektor w*, co ilustruje rysunek 


'por. Osowski [ 2], str. 26 








7.2. Zasada Hebba 
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7.2. Jeżeli wektor nie jest jednostkowy, to wartość y-, jest proporcjonalna do wielkości 
tego rzutu, albowiem: 

wfx = Yw ijXj = ||w, ; || • j|x|| • cos a = |||| • yi, 

3 S V ^ 


skąd 

E j W ij X j 


7.2.2. Sformułowanie zasady Hebba 

W neuro-fizjologii zauważono, że waga powiązań między dwoma neuronami wzrasta 
przy jednoczesnym pobudzeniu obu neuronów, w przeciwnym przypadku maleje. 

w i(k + 1) = w i(k) + Aw i(k), 

gdzie ogólnie 

A w i(k) = F{x(k),yi{k)). 

Funkcja F jest funkcją sygnału wejściowego x(/c), nazywanego również sygnałem pre- 
synaptycznym, a yi oznacza sygnał wyjściowy (pojawiający się na i-tym Wyjściu) 
nazywany również sygnałem post-synaptycznym. 

W klasycznym ujęciu Hebba funkcja F upraszcza się do funkcji iloczynowej. 
Przy uczeniu bez nauczyciela funkcja F przyjmuje wtedy postać: 


Awi(fc) = ri ■ x(fc) • yi(k ), (7.1) 

natomiast przy uczeniu z nauczycielem w powyższym wzorze zamiast sygnału post- 
synaptycznego yi{k) należy podstawić wartość docelową (target) ti{k). Formuła (7.1) 
nosi nazwę prostej reguły Hebba. 

7.2.3. Reguła Oji 

Przy użyciu funkcji iloczynowej użytej w prostej regule Hebba wagi wzrastają wy¬ 
kładniczo z postępem uczenia i wielokrotną prezentacją tego samego wzorca, co jest 
zjawiskiem bardzo niepożądanym. 

Aby temu zapobiec, wprowadzono reguły ograniczające przyrosty wag. Jedną ta¬ 
ką modyfikacją jest reguła zapominania. Inną ważną modyfikację wprowadził Oja, 
który zaproponował, aby regułę Hebba stosować na sygnale pre-synaptycznym sko¬ 
rygowanym wg następującego wzoru: 


Aw i(k) = r] ■ yi(k)[x(k ) - 2/»(fc)wi(ft)]- 


(7.2) 
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Tutaj x jest zmniejszane (proporcjonalnie) o wielkość rzutu wektora x na wektor w,;. 
Gdyby wektor w, byl jednostkowy, to byłby to dokładnie rzut (porównaj rysunek 
7.2 i uwagi na ten temat w sekcji 7.2.1). 

Oj a pokazał, że przy uczeniu określonym wzorem (7.2) wektor w,; w trakcie ucze¬ 
nia stabilizuje się i przyjmuje wartość jednostkową. 

Po wielokrotnej (kilka set lub kilka tysięcy) prezentacji próbek danych wartości 
średnie lub też wartości oczekiwane przyrostów wag powinny wynosić zero. Oznacza 
to, że po ustaleniu punktu równowagi 

£{Aw i(k)} = 0, 


i tym samym również 


' yi(k)[x(k) - yi(k)wi(k)]} = 0. 


Opuszczając współczynnik ry oraz wskaźniki (fc) pokazujące zależności poszcze¬ 
gólnych wielkości od numeru iteracji k, otrzymujemy: 


£{xyi} - £{yiViWi} = 0 . 

Biorąc pod uwagę, że yi = wfx = x T Wj, otrzymujemy w dalszym ciągu: 

£{xX T W,;} — £{wf XX^WWj} = 0. 

Oznaczmy 

£{xx t } = C, oraz wfCw, ; = \ . (7.3) 

W zależności od przeskalowania danych macierz C jest nazywana macierzą iloczy¬ 
nów skalarnych, macierzą kowariancji lub macierzą korelacji (danych). Po wykonaniu 
podstawienia określającego macierz C i skalar A,; otrzymujemy: 

Cw, = AjWj. (7.4) 

Jak widać, wektor Wj w stanie równowagi przedstawia sobą wektor własny ma¬ 
cierzy C związany z wartością własną A tej macierzy. Można pokazać, że musi 
to być wektor o długości 1: mnożąc (7.4) lewostronnie przez wf i podstawiając 
wfw i = || Wi || otrzymujemy A i = A* 11 w. ź 11 , co jest prawdziwe tylko wtedy, gdy 

ll w i|| = 1- 

Z innego rozumowania (nieprzedstawionego tutaj) wynika, że wektor w, wskazu¬ 
je na tak zwany kierunek główny zbioru danych, czyli oś główną elipsy (elipsoidy) 
przedstawiającą koncentrację punktów w przestrzeni danych. 

Oja, a później Sanger, podali efektywne wzory, jak znajdywać jednocześnie wek¬ 
tory wi,..., Wfc przedstawiające tzw., kierunki główne danych. Technika ta ma duże 
zastosowanie przy redukcji liczby cech w danych: wektory x = [x\,..., Xd\ T o dużej 





7.3. Uczenie z konkurencją 
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liczbie wymiarów można wtedy zastąpić wektorami z = [z ±,..., Zk] T , (k < d) co na 
ogół znacznie mniejszej liczbie cech. 

Osowski (str. 235) pokazuje przykład zastosowania tej techniki do zapamiętywa¬ 
nia i odtwarzania obrazów; pokazuje rn.in. obraz odtworzony na podstawie dwóch 
pierwszych składowych głównych. 


7.3. Uczenie z konkurencją 

Neurony otrzymują na początku przypadkowe wartości wag, a następnie, na podsta¬ 
wie prezentowanych danych „uczą się” rozpoznawać te dane i zbliżają się odpowiednio 
do obszarów zajmowanych przez te dane. 

Po każdej prezentacji wzorca x(fc) zwycięzcą zostaje tylko jeden neuron, najbliższy 
prezentowanemu wzorcowi (bliskość jest liczona w sensie ustalonej metryki). 

Przypuśćmy, że w tej konkurencji zwyciężył neuron o numerze c. Ma on prawo 
uaktualnić swoje wagi według jednej z dwóch zasad: WTM - Winner Takes Most, 
i WTA - Winner Takes Ali. 

• Jeżeli obowiązuje zasada WTA 1 to tylko neuron zwycięski (oznaczany dalej 
indeksem c, od conąueror ) uaktualnia swe wagi, tzn. zbliża się do wektora x(k). 
Wagi zwycięskiego neuronu mogą zostać uaktualnione np. na podstawie wzoru 
Kohonena: 

w c {k + 1) = w c (fc) + r)(k)[x(k) - w c (fc)]. (7.5) 

Współczynnik uczenia r?(fc) jest na ogól malejącą funkcją numeru iteracji k. 

• Jeżeli obowiązuje zasada WTM to również neurony sąsiadujące z neuronem 
zwycięskim, czyli należące do sąsiedztwa J\f c (k), mogą to uczynić według zasady: 

Wi(k + 1) = w.j(fc) + r/(k)G(i,c)[x(k) - w,(A:)], i € Af c (k). (7.6) 

Funkcja G oznacza tu wpływ sąsiedztwa i jest omawiana bliżej w rozdziale 8. 

W efekcie takiego współzawodnictwa następuje samoorganizacja procesu uczenia. 
Neurony dopasowują swoje wagi w ten sposób, że przy prezentacji grup wektorów 
wejściowych zbliżonych do siebie zwycięża zawsze ten sam neuron. 

Neuron, poprzez zwycięstwo we współzawodnictwie rozpoznaje swoją kategorię. 

Ogólnie można powiedzieć, że: Przy podaniu na wejście sieci wielu wektorów zbli¬ 
żonych do siebie będzie zwyciężać ciągle ten sam neuron, w wyniku czego jego wagi 
będą odpowiadać uśrednionym wartościom wektorów wejściowych, dla których dany 
neuron byl zwycięzcą. 

Neurony nie wy grywające nie zmieniają swoich wag. Mówi się, że pozostają 
one martwe. 
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7.3.1. Przykład uczenia WTA wektorów leżących na kole 

Zasadę uczenia z konkurencją WTA ilustruje rysunek 7.3 zaczerpnięty z książki 
Osowskiego. Zbiór danych składa się z ośmiu wektorów; każdy z tych wektorów jest 
dwuwymiarowy i ma długość jednostkową, wobec tego przestrzeń danych może być 
odwzorowana na okręgu (w przypadku wektorów danych wielowymiarowych o dłu¬ 
gości 1 - wektory danych mogłyby być odwzorowane na sferze). 

Mamy również 4 neurony z przypisanymi do nich wektorami wagowymi W; = 
[wii,Wi 2 ) T , i = 1,2,3,4, każdy z nich o długości jeden, tj. ||w*|| = 1, Vj. 

Neurony te miały się, w procesie konkurencyjnym uczenia, nauczyć prezentowa¬ 
nego zbioru danych i zaadaptować swoje wagi w ten sposób, aby możliwie dobrze 
reprezentować te dane. 

Uczenie było konkurencyjne (WTA), a elementy zbioru danych prezentowano neu¬ 
ronom w sumie 320 razy. 



Rysunek 7.3: Uczenie się z konkurencją. Jeden neuron nigdy nie zwyciężył , wskutek 
czego nie zdołał wytworzyć swojej strefy wpływów 


Po prezentacji kolejnego wzorca zwycięzcą zostawał neuron najbliższy temu wzor¬ 
cowi (najbliższy w sensie odległości kątowej mierzonej po okręgu koła). 

Po zakończeniu prezentacji okazało się, że w poszczególnych prezentacjach zwy¬ 
cięzcami zostawały tylko trzy neurony; i tylko one uzyskały przywilej zaadaptowania 
swoich wag. Jeden z neuronów nigdy nie zwyciężył, wobec czego jego wagi pozo¬ 
stały niezmienione. Pozostałe neurony uzyskały swoje strefy wpływów, W rezultacie 
neurony te uplasowały się w okolicy średnich pewnych grupek danych i stały się 
reprezentantami tych grupek danych. 









7.3. Uczenie z konkurencją 
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7.3.2. Kwantyzacja przestrzeni danych i wieloboki Voronoia 

Opisany w poprzedniej sekcji przykład uczenia się neuronów dla danych należących 
do R 2 może być zrealizowany dla danych należących do dowolnego zbioru punk¬ 
tów w R d . W trakcie procesu uczenia neurony starają się uplasować pośród tych 
punktów-danych i stać się ich reprezentantami. W ten sposób cala przestrzeń zosta¬ 
je podzielona na obszary atrakcji odzwierciedlające strefy wpływów poszczególnych 
neuronów. W przypadku używania metryki eukliclesowej są to obszary wypukłe, które 
są nazywane są czasem wielobokami Voronoia, a reprezentujące je wektory — wek¬ 
torami Voronoia. Obszary Voronoia Vi są zdefiniowane jako miejsce geometryczne 
punktów x spełniających następujący warunek: 

Vi = {x: || w* - x1 1 < ||wj - x||, j} 

Przykład wieloboków Voronoia jest pokazany na rysunku 7.4 — gdzie mamy prze¬ 
strzeń danych określoną jako prostokąt. 



Rysunek 7.4: Wieloboki Voronoia utworzone przez 7 wektorów kodowych na prosto¬ 
kącie z metryką euklidesową 


Podkreślmy tutaj, że używanie różnych metryk (tj. korzystanie z różnych definicji 
odległości między dwoma punktami leżącymi w przestrzeni danych) daje różne po¬ 
staci obszarów atrakcji. Przykład rozkładu takich obszarów atrakcji — przy użyciu 
czterech różnych metryk — jest pokazany w książce Osowskiego, str. 252. 

Kohonen nazwał proces tworzenia reprezentantów danych kwantowaniem wek¬ 
torowym ( Vector Quantization) , lub dokładniej: adaptacyjnym kwantowaniem wek¬ 
torowym ( LVQ, Learning Vector Quantization). Wektory wagowe neuronów zostały 
przez Kohonena nazwane słowami kodowymi, codebook nectors , a ich zbiór — książką 
kodową, codebook 

W ten sposób cała przestrzeń danych została skwantowana; a wektory danych 
należące do jednego wieloboku Voronoia można zastąpić odpowiadającym mu słowem 
kodowym. 







72 
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Osowski (str. 268-275) podaje następujące przykłady zastosowań sieci samoorga¬ 
nizujących się: 

1. Kompresja obrazów 

2. Wykrywanie typu uszkodzeń 

3. Prognozowanie obciążeń systemu elektroenergetycznego 



Rozdział 8 


SOMy czyli samoorganizujące 
się mapy 


Typowy przedstawiciel takich sieci: Sieci Kohonena 1 . 

SOM-y (Self Organizing Maps) realizują generalnie dwa zadania: 

1. Wektorowej kwantyzacji (kompresji danych), 

2. Odtwarzanie przestrzennej organizacji danych wejściowych 

8.1. Oznaczenia 

Niech x = [xi,... ,Xd] T oznacza d -wymiarowy wektor danych, tzw. próbkę lub wzo¬ 
rzec. Wektor x może być interpretowany jako punkt d-wymiarowej przestrzeni: x € 
R d . 

Załóżmy, że mamy to neuronów tworzących mapę. Neurony te mogą być ułożo¬ 
ne w siatkę (ml x m2) na płaszczyźnie, lub liniowo na odcinku 1 x to). Możliwe 
(ale znacznie rzadziej wykorzystywane) są siatki trójwymiarowe, lub rozłożone na 
cylindrze lub torusie. 

Każdy neuron jest scharakteryzowany swoim wektorem wag (nazywanym również 
wektorem kodowym lub wektorem Voronoia). Tak więc mamy m wektorów wagowych 
wi,... ,w m , gdzie w* = [w ily Wi 2 ,. . - ,w id \ T , i = 1, ■. •, to, w* G R d . 

Te wszystkie pojęcia występowały już wcześniej. 

Teraz definiujemy dodatkowe wektory referencyjne ri,.. .,r m określające po¬ 
zycje węzłów siatki. Wektory referencyjne r*, (i = 1,..., to), są ściśle przyporządko¬ 
wane wektorom kodowym Wj należącym do przestrzeni danych wejściowych (input 
space) R d . 

1 opisane np. w książce Osowskiego, str. 249-275 



74 


Rozdział 8. SOMy czyli samoorganizujące się mapy 


W dalszym ciągu będziemy zajmować się tylko mapami na płaszczyźnie, wobec 
tego r, e R 2 , i = 1,..., m. 

Inicjacja wektorów kodowych i przyporządkowań do nich odpowiednich wektorów 
referencyjnych jest na ogól przypadkowa. Sieć (mapa) na podstawie prezentowanych 
wzorców powinna nauczyć się danych. Uczenie jest konkurencyjne. 

Jak wiemy, przy uczeniu konkurencyjnym zwycięża tylko jeden neuron (jego nu¬ 
mer oznaczaliśmy indeksem c), ale neurony znajdujące się w sąsiedztwie neuronu- 
zwycięzcy mogą mogą też częściowo partycypować w zwycięstwie neuronu c i adap¬ 
tować swoje wagi. Dlatego niezmiernie ważnym pojęciem jest pojęcie sąsiedztwa 
zwycięskiego neuronu. 


8.2. Siatki i sąsiedztwa 


Istotą sieci Kohonena jest ustalenie pewnej sztywnej siatki neuronów. Siatka taka 
(grid, map) ma najczęściej budowę hexagonalną lub prostokątną, jak to pokazano na 
rysunku 8.1. 

Na siatce heksagonalnej każdy neuron ma sześć, a na prostokątnej osiem sąsiadów 
pierwszego rzędu. Można by liczyć również sąsiedztwo 2-go rzędu, i dalsze, jak to 
pokazano na rysunku 8.1. 
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(b) Rectangular grid 


Rysunek 8.1: Siatka hexagonalna i prostokątna 


Sąsiedztwo można wyznaczać również analitycznie, za pomocą pewnych funkcji 
wyznaczających odległość neuronu i od zwycięskiego neuronu c. Funkcje te mają na 
ogól postać funkcji radialnych, zcentrowanych w wektorze referencyjnym r c . 

Najbardziej popularnymi funkcjami sąsiedztwa są bubble i gaussian. Są one poka¬ 
zane na rysunku 8.2. Przykładowo funkcja gaussian przyjmuje postać 

hgau SS (i,c) = exp 

Jeszcze inne funkcje sąsiedztwa ( cut-gaussian i Epanechnikov ) są pokazane na rysun¬ 
ku 8.3. 

Sąsiedztwo neuronu c będziemy oznaczać A f c (od neighbourhood). Jeżeli chcemy 
wyraźnie napisać, że jest to sąsiedztwo neuronu c który zwyciężył w iteracji k, to 


( Ilrc-Ulh 
l 2 * 2 (t) ) 

















8.3. Adaptacja wektorów kodowych podczas procesu uczenia 
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Rysunek 8.2: Funkcje sąsiedztwa: bubble i gaussian 



Rysunek 8.3: Jednowymiarowe funkcje sąsiedztwa dla promienia R=2: bubble, gaus¬ 
sian, cut-gaussian, Epanechnikov 


zapiszemy 

A f c {k). 

Na ogół sąsiedztwo Af c (k) charakteryzuje się pewnym promieniem, który maleje 
wraz z upływem czasu uczenia t, czyli w miarę zwiększania się wskaźnika k, oznacza¬ 
jącego numer kolejny iteracji (prezentacji wzorców próbki uczącej). 


8.3. Adaptacja wektorów kodowych podczas procesu 
uczenia 

Proces uczenia polega na prezentowaniu sieci wzorców zbioru uczącego, z generalnym 
celem, żeby siec upodobniła swoje wektory kodowe do prezentowanych wzorców. 
Podczas uczenia obowiązuje konkurencja. 

Niech k będzie wskaźnikiem (numerem) kolejnej prezentacji. Po prezentacji sieci 
wektora danych x(fc) zwycięża ten neuron, którego wektor kodowy jest najbliższy 
dopiero co zaprezentowanemu wektorowi x(fc) w sensie ustalonej metryki (najczęściej 
jest to metryka euklidesowa). 

Oznaczmy numer neuronu-zwycięzcy symbolem c (conąueror). Wektor kodowy 
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tego neuronu spełnia relację: 


d(x,w c ) = min d(x, Wj). 

Zwycięski neuron oraz neurony należące do jego sąsiedztwa (o ile aktualny zbiór 
Af c {k) obejmuje jakichś sąsiadów) mogą aktualizować swoje wagi, to znaczy upodab¬ 
niać je do aktualnego zaprezentowanego wzorca x(fc), natomiast wagi neuronów nie 
należących do A f c (k) nie zmieniają się. 

8.4. Formuły na zmianę wag wygrywających neuro¬ 
nów 

Jeżeli i-ty neuron otrzymał prawo zmiany swoich wag, to zmiana ta dokonuje się 
według wzoru: 


w i(k + 1) = w i(k) + r](k) ■ hd(k) • [x(fc) - w,(fc)], 


gdzie f)(k) oznacza współczynnik uczenia — piszemy o nim w sekcji 8 . 6 . 
c oznacza numer wektora-zwycięzcy, tj. znajdującego się najbliżej prezentowanego 
w fc-tym kroku wektora x(fc), 

h c i(k) określają, czy neuron i należy do sąsiedztwa zwycięskiego neuronu. 

Jeśli chodzi o współczynnik uczenia rji(k), to na ogół maleje on wraz z upływem 
czasu uczenia wyznaczanego numerem iteracji k (por. sekcja 8 . 6 ); 

Jeśli chodzi o funkcję h C i(k ), to może ona być określona w następujący sposób: 


h C i{k) — 


1 gdy i € N c (k ), 
0 gdy i i N c (k). 


Funkcja h C i(k) może też zależeć bezpośrednio od odległości D między wektorami 
referencyjnymi r c i rą, np. 


h c i{k) 


g(D(c,i)) gdy i G N c (k), 
0 gdy i i. N c (k ), 


gdzie g(-) jest funkcją malejącą (dokładniej: nierosnącą) przy wzroście swego argu¬ 
mentu. Może to być np. funkcja bubble lub gaussian pokazana na rysunku 8.2. 


8.4.1. Dwie fazy uczenia 

Na ogół uczenie przebiega w dwóch fazach. Najpierw przyjmuje się duże wartości 77 
i duży promień sąsiedztwa. 

W drugiej fazie (fine tuning) obydwie te wielkości ulegają istotnemu zmniejszeniu; 
w szczególności promień sąsiedztwa maleje do zera. 




8.5. Uczenie sieci on-line i wsadowo, współczynnik uczenia 
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Pierwsza faza - przebiega według zasady WTM - promień sąsiedztwa jest duży, 
co powoduje, że oprócz neuronu-zwycięzcy również jego sąsiedzi z mapy zmieniają 
swoje wektory kodowe. Również współczynnik uczenia r/ jest stosunkowo duży. 

Tak więc, przy każdej prezentacji kolejnego wektora x zostanie do niego przycią¬ 
gnięty odpowiadający mu wektor-zwycięzca, który pociąga za sobą neurony z naj¬ 
bliższego sąsiedztwa na zdefiniowanej siatce. 

Jeśli prezentujemy sieci np. wektory xi, X 2 , X 3 położone blisko siebie w przestrzeni 
R N , to znajdą one ten sam neuron w który będzie pociągał za sobą swoje otoczenie. 

Druga faza uczenia. Obowiązuje tu zasada WTA. Adaptacji podlegają tylko 
neuron-zwycięzca c, ponieważ promień sąsiedztwa zmalał do zera. Wzór na adapta¬ 
cję: 

w c (k + 1) = w c (k) + 77(fc)[x - w c (fc)]. 

Ostatecznie cała przestrzeń R d w której znajdują się dane, zostanie podzielona 
na strefy wpływów poszczególnych neuronów (obszary Voronoia). Sąsiedztwo wektory 
referencyjnych na mapie powinno odzwierciedlać (topologicznie) sąsiedztwo wektorów 
kodowych w przestrzeni danych. 


8.5. Uczenie sieci on-line i wsadowo, współczynnik 
uczenia 

Uczenie może się dokonywać na dwa sposoby: sekwencyjnie i wsadowo (batch). 

Uczenie sekwencyjne , inaczej na bieżąco, lub on-line polega na tym, że dla t = 
1,2,... prezentujemy sieci wektory danych x(t), po czym następuje uaktualnienie 
wag zwycięskiego neuronu (i ewentualnie jego sąsiadów), jak to pokazano w sekcji 
8.4. 

Współczynnik uczenia na ogól zmniejsza się z czasem. 

Niech T oznacza maksymalną liczbę iteracji. Stosuje się następujące wzory na 
zmniejszanie współczynnika uczenia: 

1. Liniowe zmniejszanie 

= t = 1,2,.... 

2. Wykładnicze zmniejszanie 

r](t) = rj 0 exp(—Ct), t = 1,2,..., C > 0 jest pewną stalą. 

3. Hiperboliczne zmniejszanie 

rj(t) = , t = 1,2,..., Ci, C% > 0 pewne stałe. 

4. Indywidualny współczynnik uczenia, np. 

r]i(t) = l/rii(t), gdzie rij(f) oznacza liczbę zwycięstw i-tego neuronu. 

Wsadowe uczenie w przypadku SOM-ów jest znacznie szybsze i bardziej sta¬ 
bilne; polega na wykonywaniu aktualizacji wag na podstawie wszystkich próbek. Jest 
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to wariant domyślny w pakiecie somtoolbox vs.2. Algorytm uczenia wsadowego jest 
następujący (por. Skubalska, [1 ], za Kohonenem [7]) 

1. Ustal początkowe wektory kodowe i odpowiadające mu wektory referencyjne na 
mapie. Takimi początkowymi wektorami kodowymi może być np. m dowolnych 
wektorów danych z próbki uczącej. 

2. Dla każdego wektora w, zapamiętaj zbiór wektorów uczących {x^} które 
oddziaływałyby na w, w zwykłym algorytmie SOM. 

3. Wyznacz nowe wartości wag (symbol c(xj,) oznacza neuron wygrywający przy 
prezentacji wektora danych x^, natomiast N jest ogólną liczebnością próbek 
uczących) 


N 


N 



4. Jeśli nie jest spełnione kryterium STOP-u (tutaj nie podaliśmy go), wróć do 
kroku 2. 


8.6. Dobroć aproksymacji 

Błąd, kwantyzacji — jako średnia odległość punktów od ich reprezentantów czyli wek¬ 
torów wag (data representation accuracy, average ąuantization error between 
data vectors and their BMU - best maching units — on the map). 

Błąd topologicznej reprezentacji — procent wektorów danych, dla których pierwsi dwaj 
najbliżsi reprezentanci nie są sąsiadami na siatce (data set representation ac¬ 
curacy, the topographic error; percentage of data vectors for which the first- 
and second-BMU are not adjacent units). 

8.7. Dostępne dla nas oprogramowanie 

SOM_PAK zestaw programów w języku C autorstwa Kohonena i jego zespołu. 
Jest dostępny w internecie, tworzy grafikę w postaci .ps lub .eps; wymaga kom¬ 
pilacji (make). Pracuje zarówno na Unix-ie jak i na PC-tach. Adres internetowy: 
http://www.cis.hut.fi/research/som_lvq_pak 

somtoolbox vs. 2 ([18]) - zestaw skryptów w postaci M-files autorstwa J. Vesanto 
i współautorów. Dostępny w internecie pod adresem : 
http://www.cis.hut.fi/projects/somtoolbox 






8.8. Mapa 49 województw polskich 


wizualizacja U-mat 
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8.8. Przykład tworzenia mapy dla 49 województw 
polskich; wizualizacja U-mat 

Rozpatrujemy tablicę danych o wymiarze 49 x 9, której wiersze odpowiadają 49 
województwom polskim (dane pochodzą z r. 1990), z których każde zostało scharak¬ 
teryzowane przez 9 cech socjo-ekonomicznych (kolumny tablicy). Tym samym każde 
województwo może być interpretowane jako punkt w R 9 . 

Dla grupowania punktów zastosowano mapę heksagonalną o siatce 6x6. 

Na rysunku 8.4 pokazujemy mapę otrzymaną w wyniku obliczeń oryginalnym 
programem SOMLPAK Kohonena (niestety, program ten nie liczy błędu topologicznej 
reprezentacji). 



Rysunek 8.4: Mapa Kohonena obrazująca zgrupowania 49 województw polskich ze 
względu na 9 cech socjo-ekonometrycznych. Niektóre węzły mapy pozostały puste, 
inne wektory kodowe zdołały przyciągnąć po kilka punktów-województw 


To co widzimy na siatce (mapie) jest obrazem wielowymiarowej przestrzeni. Siat¬ 
ka składa się z obszarów heksagonalnych, w których środkach znajdują się wektory 
referencyjne odpowiadające wektorom kodowym umiejscowionym w R 9 . Oznacza to, 
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że każdemu neuronowi r,; na mapie odpowiada wektor kodowy (codebook vector) w,; 
w przestrzeni wejściowej (input space) R 9 . 

Faktyczne odległości między wektorami w* są obrazowane odcieniami szarości na 
mapie: obszary bliskie są jasne, ciemny kolor oznacza duże odległości, a więc może 
oznaczać granice klasterów z R 9 (inne programy, np. Somtoolbox operują w tym 
celu kolorami z odpowiednim kluczem na oznaczenie bliskich i dalekich wektorów 
wagowych). Na utworzonej mapie właściwe heksagony (j) zawierające węzły mapy są 
otoczone dodatkowymi heksagonalni pokazującymi kolorystycznie, jaka jest średnia 
odległość wektora kodowego (j) od sąsiadujących z nim wektorów kodowych. 

Taki sposób wizualizacji nosi nazwę u-mat od informatyka o nazwisku Ultsch, 
który zaproponował taki sposób wizualizacji map. 

Mając klucz do województw możemy próbować interpretować powstałe zgrupo¬ 
wania. Punkty 1 i 24 to województwa Warszawskie i Łódzkie. Punkty 18 i 47 to 
Kraków i Wrocław; punkt 32: Poznań; punkt 3 i 22 to Białystok i Lublin; punkt 10: 
Gdynia-Gdańsk. 

Wszystkie te punkty to miasta uniwersyteckie z pewną tradycją. Tworzą one 
wyraźny klaster, oddzielony od pozostałych punktów zatoką pustych węzłów. 

Innym widocznym na mapie klasterom można przypisać również ciekawą inter¬ 
pretację. 


8.9. Inne algorytmy uczenia sieci samoorganizują¬ 
cych 

Miary odległości między wektorami, problem normalizacji wektorów (należy dodać 
tzw. normalizację na 0-1, o czym Osowski nie mówi) 

• Algorytmy uwzględniające „zmęczenie” neuronów po zwycięstwie. 

• Algorytm stochastycznej relaksacji (neurony ulegają adaptacji z pb-stwem okre¬ 
ślonym wzorem Gibbsa). 

• SCS, Soft Competition (uwzględnia się stopień aktywności neuronów, fawory¬ 
zując jednostki mało aktywne). 

Współczynnik uczenia może być indywidualny dla każdego neuronu, jest uza¬ 
leżniony od tego, co robił neuron w ostatnich k — 1 prezentacjach 

• Algorytm gazu neuronowego. Podobno bardziej skuteczny niż algorytm Koho- 
nena, który jest jego szczególnym przypadkiem. 

Neurony są sortowane w zależności od ich odległości od wektora x a następ¬ 
nie ulegają adaptacji w zależności od ich pozycji s(?') w posortowanym ciągu 
odległości - wg funkcji wykładniczej exp(—s(i)/ A). 



Rozdział 9 


Mapy SOM — Pakiet 
SomToolbox2 


Pakiet ten działa w środowisku MATLAB 5. Podstawą obliczeń są dwie struktury: 
data-struct (sD) i som-strukt (sM). 

9.1. Tworzenie struktury danych data-struct sD i nor¬ 
malizacja 

Strukturę taką można utworzyć na dwa sposoby: 

1. sD = som_read_data(’iris.dat'); - przez wczytanie danych z pliku Ascii; 
funkcja może mieć dodatkowy argument (,’x’) określający braki w danych czyli 
missing values , 

2. sD = som_data_struct(D, 'namej ’iris-sd’, 'comp-names', {'SepalL, ... ,'PetalW'} ) 

; - przez utworzenie struktury sD z MatLabowskiej tablicy danych znajdującej 
się w przestrzeni roboczej (workspace) MatLaba. 

O polach utworzonej w ten sposób struktury sD piszemy dalej w sekcji 9.1.3. 

9.1.1. Pierwszy sposób — czytanie danych Ascii 

Dane w pliku (np. o nazwie ’iris.data’) powinny mieć następującą postać (kropki 
oznaczają, że opuściliśmy tu pewne fragmenty tych danych): 

4 

#n SepalL SepalW PetalL PetalW 
5.1 3.5 1.4 0.2 Setosa 
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4.9 3.0 1.4 0.2 Setosa 

5.0 3.3 1.4 0.2 Setosa 
7.0 3.2 4.7 1.4 Versicolor 


5.9 3.0 5.1 1.8 Virginica 

Strukturę Data-struct zawierającą te dane można utworzyć wykonując instrukcję: 
sD=som_read_data('iris.data'); 

Gdybyśmy nie mieli w pliku nazw osobników (nazw przypisanych poszczegól¬ 
nym wektorom-wierszom tablicy danych), to możemy je dodać później do utworzonej 
struktury sD za pomocą instrukcji som Jabel - opisanej w następnej sekcji. 

9.1.2. Drugi sposób — korzystanie z tablicy MatLabowskiej D 

Jeśli mamy tablicę D znajdującą się w przestrzeni roboczej MatLaba, to możemy 
strukturę sD utworzyć kolejno za pomocą instrukcji: 

sD = som_data_struct(D, , name’ 1 ’iris-sD l , 'compmames', 'SepalL’, 'SepalW', ’PetalL’,’PetalW') 

Utworzona w ten sposób struktura danych zawiera nazwy zmiennych (’SepalL, ... 

, ’PetalW’), ale nie zawiera nazw osobników, tj. wektorów wierszy. Można je dodać 
za pomocą rozkazu somJabel (tutaj każdy wektor-wiersz otrzymuje po prostu nazwę 
gatunku irysa do którego należy): 

sD = somJabel(sD,'add’,[1:50]’,'Setosa'); 
sD = somJabel(sD,'add’, [51:100]', ’Versicolor'); 
sD = som _la bel (sD, 'add’,[101:150]’,’Virgi n ica ’); 

9.1.3. Pola struktury danych sD 

Utworzona instrukcjami przedstawionymi w sekcji 9.1.1 lub 9.1.2 struktura sD zawiera 
następujące pola: 

sD = 

type: ’som_data’ 
data: [150x4 double] 
labels: {150x1 celi} 
name: ’ iris.data’ 
comp_names: {4x1 celi} 
comp_norm: {4x1 celi} 
label_names: [] 

Pole comp_names zawiera nazwy zmiennych. 

Pole compmorm zawiera informacje o normalizacji danych. 




9.1. Tworzenie struktury danych data-struct sD i normalizacja 
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9.1.4. Normalizacja i denormalizacja zmiennych 

Dopuszcza się następujące możliwości normalizacji: 

’var’ (na 0-1), tj. odjąć średnią i podzielić przez a, 

’range’, tj. na min-max, 

’log’ ( x new = ln(x — min(x) + 1)), 

logistic’ (softmax; x = (x — x)/a x \ x tr = 1/(1 + exp{—x}), 

IristD’, ’histC’ (histogram eąualization). 

Normalizacja danych wykonuje się za pomocą instrukcji: 
sD = som_normalize(sD,'var’); % zamiast ’var’ może być inny sposób, np. ’range’, 
log’, ... . 

Funkcja sommornralize może zawierać jeszcze trzeci argument określający które 
zmienne mają być normalizowane. Np. 

sD = som_normalize(sD,’log',[l 3]); - zostaną zlogarytmowane tylko pierwsza i trzecia 
zmienna. 

Przykładowa informacja zawarta w 1 komórce pola sD.compmorm: 


sM.comp_norm{l} 
ans = 

type: ’som_norm’ 
method: ’var’ 
params: [5.8433 0.8281] 
status: ’done’ 


Normalizacja wykonana na strukturze sD może zostać wykonana w ten sam sposób 
na innej (nowej) macierzy o nazwie Dn. Należy wydać w tym celu rozkaz: 

Dn = som_normalize(Dn, sD); 


Denormalizacja danych jest wykonywana za pomocą analogicznych instrukcji: 
sD = som_denormalize(sD); 

Zostaną unieważnione wszystkie normalizacje wykazane w polu .compmorm i przy¬ 
wrócone wartości danych sprzed normalizacji. 

Zwróćmy uwagę, że dane zostają „zdenormalizowane”, ale opis normalizacji po¬ 
zostaje. Aby usunąć opis normalizacji należy użyć przy denormalizacji kwalifikatora 
’remove’: 

sD = som_denormalize(sD,'remove'); 
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9.2. Tworzenie mapy, instrukcja som_make 

9.2.1. Postępowanie standardowe — instrukcja som_make z war¬ 
tościami domyślnymi 

Struktura-mapa sM może być utworzona za pomocą rozkazu (zaleca się, żeby dane 
zostały najpierw znormalizowane): 

sM = som_make(sD); 

Procedura som_make wywołana bez dalszych parametrów inicjalizuje i trenuje 
mapę według wartości domyślnych wbudowanych w procedurę. 

W ten sposób zostają automatycznie określone rozmiary mapy (na podstawie 
stosunku wartości własnych m. kowariancji obliczanych danych) a następnie (default) 
wykonuje uczenie wsadowe (training using batch algorithm) w dwóch fazach (Rough 
training phase... Fine tuning phase...), po czym zostaje obliczony błąd kwantyzacji 
i błąd reprezentacji topologicznej (opisane wcześniej w sekcji 8.6). 

Te ostatnie mogą zostać obliczone oddzielnie za pomocą rozkazu: 

[q,t] = som_quality(sM,D) % for the trained map sM 
Dla danych Iris otrzymujemy: 

Finał ąuantization error: 0.393 Finał topographic error: 0.013 


9.2.2. Pola struktury map-struct 

Utworzona struktura sM zawiera następujące pola: 


sM = 

type: 
codebook: 
topol: 
labels: 
neigh: 
mask: 
trainhist: 
name: 
comp_names: 
comp_norm: 


’som_map’ 

[66x4 double] 

[1x1 struct] 

{66x1 celi} 
'gaussian’ 

[4x1 double] 

[1x3 struct] 

’SOM 29-Dec-2000’ 
{4x1 celi} 

{4x1 celi} 


Pole codebook zawiera wektory wag; pole labels jest na razie puste (wypełni¬ 
my je za chwilę, wykonując som_autolabel), comp_names zawiera nazwy zmiennych, 
comp_norm - informacje i dane do normalizacji. 

Pole codebook zawiera współrzędne wektorów wagowych: 

sM.codebook 
ans = 

-1.4152 0.0229 -1.3285 -1.3285 
-1.3005 0.0068 -1.2852 -1.2746 




9.2. Tworzenie mapy, instrukcja som make 
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1.2796 0.4892 1.1476 1.3781 

1.6463 0.5612 1.2880 1.3040 

Pole topol zawiera informacje o topologii utworzonej mapy: 

sM.topol 
ans = 

type: ’som_topol’ 
msize: [11 6] 
lattice: ’hexa’ 
shape: ’sheet ’ 

Pole labels może zawierać nazwy wektorów wagowych. Na początku, po utworze¬ 
niu mapy, pole to jest puste. Odpowiednie nazwy można nadać za pomocą procedury 
som_autolabel patrz niżej. 

Pole mask zawiera tzw. maskę na zmienne: jeśli elementem maski jest wartość 
zerowa, to odpowiednia zmienna nie jest uwzględniana w obliczeniach. 

sM.mask’ 

ans =1111 

compnames to oczywiście nazwy zmiennych: 
sM.comp_names’ 

ans = 'SepalL’ 'SepalW’ ’PetalL’ 'PetalW’ 

Pole comp.norm zawiera informacje o normalizacji, takie same jak struktura sD 
omawiana wcześniej (por. sekcja 9.1.3). 

9.2.3. Postępowanie niestandardowe — instrukcja som_make 
z deklarowanymi wartościami parametrów 

Ogólna postać wywołania funkcji sormmake: 
function sMap = som_make(D, varargin); 

Dla możliwych argumentów najpierw podaje się identyfikator argumentu, a potem 
jego wartość. Przykłady: 

sMap = som_make(D, [[argID,] value, ...]); 

sMap = som_make(D); 

sMap = som_make(D, 'munits', 20); 

sMap = som_make(D, 'munits', 20, ’hexa', 'sheet'); 

sMap = som_make(D, 'msize', [4 6 7], 'lattice', 'rect'); 

Jako D może wystąpić (matrix) training data, size dlen x dim lub (struct) data 
struct. 

Dalszymi argumentami mogą być: 
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’init’ *(string) initialization: 'randinit’ or 'lininit’ (default) 

'algorithm’ *(string) training: ’seq’ or ’batch’ (default) or 'sompak’ 

'munits’ (scalar) the preferred number of map units 

’msize’ (vector) map grid size 

'mapsize’ *(string) do you want a~’small’, 'normal’ or ’big’ map 

Any explicit settings of munits or msize override this. 
'lattice’ *(string) map lattice, ’hexa’ or ’rect’ 

’shape’ *(string) map shape, ’sheet’, ’cyl’ or ’toroid’ 

’neigh’ *(string) neighborhood function, 'gaussian’, ’cutgauss’, 

’ep’ or 'bubble’ 

’topol’ *(struct) topology struct 

’som_topol ’,’ sTopol’ = ’topol’ 

’mask’ (vector) BMU search mask, size dim x 1 

’name’ (string) map name 

’comp_names’ (string array | cellstr) component names, size dim x 1 
'tracking’ (scalar) how much to report, default = 1 

'training’ (string) ’short’, 'default’, ’long’ 

(vector) size 1x2, first length of rough training in epochs, 
and then length of finetuning in epochs 

9.3. Procedury som la bel i som autolabel 

9.3.1. Procedura somJabel 

Formalny nagłówek: som Jabel(sTo, modę, inds, [labels] ) 

Przykłady zastosowań: 

sM = som_label(sM, 'add', [1; 10], ’x’); % Dodanie próbkom nr. 1 i 10 nazwy ’x’ 
sD = som_label(sD, 'elear', 'all'); % usuwa wszystkie nazwy z danych 
sD = som_label(sD, 'replace', [1:10]', 'topten'); 

% zamienia nazwy próbek nr 1 10 na nazwę ’topten’ 

Ostatnia instrukcja mogłaby być zastąpiona następującymi dwoma: 

sD = som_label(sD, 'elear', [1:10]'); sD = som_label(sD, 'add', [1:10]’, 'topten'); 

Jeszcze inne zastosowanie: usunięcie pustych etykiet z wszystkich jednostek mapy: 

sD = som_label(sM, 'prune', 'alt); 

9.3.2. Procedura som_autolabel 

Ogólna postać procedury: som_autolabel(sTo, sFrom, [modę], [inds]); 

znaczenie: dokąd wstawić, skąd pobrać, sposób, wskaźniki 

modę: ’add’ — po prostu dodać, mogą się powtarzać 
’addl’ — zostaje zapamiętana tylko 1 etykieta 
’freq’ — dla powtarzających nazw zostaje zapamiętana 
tylko jedna nazwa i"jej frekwencja 



9.4. Wizualizacja mapy: procedura som_show 
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’vote’ — zostaje dodana nazwa o~najw. częstości 
w~przyp. losowania pierwsza wylosowania 

Wymienione operacje nie zmieniają starych nazw znajdujących się w sTo. 
Przykłady stosowania: 

sM = som_autolabel(sM,sD) % oznakowanie wg danych 

sD = som_autolabel(sD,sM) % oznakowanie wg mapy (codebook vectors) 

sD = som_autolabel(sM,sD,’vote',[5]) % etykietuje jednostki mapy 

według nazw występujących w 5-tej kolumnie danych, na zasadzie ’vote’. 

% może sie przydać, usuwa z pokazanej mapy wszystkie markery 
% (hits, labels trajectories) naniesione przez som_show_add 
% i można zaczac pokazywanie od nowa 

som_show_clear 


9.4. Wizualizacja mapy: procedura som.show 

9.4.1. Procedury som_show i som_show_add 

Standardowe wywołanie: som_show (sM); Wykreśla plansze typu U-mat i mapy w po¬ 
staci heksagonów jeżeli nie utworzono mapy o innej topologii. Default: dla wszyst¬ 
kich zmiennych. Można uzyskać wykresy tylko dla niektórych zmiennych, np. dwie 
mapy u_mat możemy otrzymać za pomocą rozkazu 
som_show(sM, 'umat', {[1 2],'1,2 only’},'umat’,{[3 4],'3,4 only’} ); 

Przykładowe wykresy są przedstawione na rysunku 9.2 a skrypt, za pomocą któ¬ 
rego wyprodukowano te wykresy, jest podany w sekcji 9.4.2. 

Można używać zestawów kolorów: (colormap) jet, hsv, hot, gray 

Metoda wizualizacji u_mat 

Metoda ta pochodzi od informatyka o nazwisku Ultsch. Na mapie widzimy podwój¬ 
ne heksagony. Struktura takiej mapy jest pokazana na rysunku 9.1. Pokazano tam 
fragment mapy hexagonalnej o wymiarach 5x3. 

Mapa U-mat pokazuje kolorami: 

a) w głównych sześciobokach oznaczonym pojedynczymi liczbami: średnie lub me¬ 
diany odległości od sąsiadów; 

b) w dodatkowych sześciobokach: odległości od wschodniego lub południowego są¬ 
siada. 

Na wykresie U-mat wysokie wartości oznaczają granice klasterów; niskie wartości 
rozłożone obok siebie wskazują na same klastery. Przykładowy wykres średnich od¬ 
ległości między sąsiadami każdego węzła mapy pokazano na rysunku 9.4 jako sub- 
plot(2,2,2). 
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Rysunek 9.1: zasady wizualizacji u_mat. 


Wykresy typowej mapy składającej się z heksagonów 

Ten typ wykresu jest pokazany na rysunku 9.2 jako wykresy górny prawy i obydwa 
dolne. 

Na tego typu wykresach można nanosić różnego rodzaju informacje o wektorach 
danych które znalazły się w strefie wpływów wektora kodowego stojącego w odpo- 
wiedniości z danym heksagonem. 

Przykładowo wykres dolny lewy na rysunku 9.2 pokazuje gatunki i frekwencje 
irysów należących do danego heksagonu, a wykres dolny prawy pokazuje kolorem 
i wielkością gatunek i frekwencję odpowiednich irysów. 

Wykres górny prawy na tym samym rysunku pokazuje tzw. component-plane dla 
pierwszej cechy irysa ( SepalL ). Tutaj dla każdego heksagonu jest pokazana średnia 
wartość tej cechy (tj. SepalL ) irysów skupionych wokół wektora kodowego będącego 
w odpowiedniości z danym heksagonem. 

Wykresy typu component-plane można otrzymać dla każdej lub dla wybranych 
zmiennych. 


Procedura som_show_add 

nanosi na wykres — wyprodukowany przez som_show dodatkowe informacje dotyczące 
etykiet (labels), obsadzenia jednostek (hits) i trajektorii (trajectories). 

labels procedura sormautolabel, nadaje etykiety jednostkom, is used to categorize 
the units (or some units) by giving them names, 
hits Hit histograms are actually markers that show the distribution of the best 
matching units for a given data set (procedura somdrits). 
trajectories - show the best matching units for a given data set that is time series 
(or any ordered series), procedura som_trajectory. 





9.4. Wizualizacja mapy: procedura som_show 
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Przykłady użycia 

som_show(sM, 'urnat', ’a11’, 'comp', 1:4, 'empty', 'Labels’, 'norm',’d'); 
som_show_add('laber,sM, 'subplot',6); 


9.4.2. Przykładowy skrypt wykorzystujący procedurę som_show 

Pokazujemy skrypt który rysuje różne warianty map Kohonena dla danych iris.data 
znajdujących się w kartotece somtb2 

•/.•/. job_soml.m 

sD=som_read_data(’iris.data’); sD=som_normalize(sD,’var’) 
sD=som_label(sD,’replace’,[1:50],’Se’); 
sD=som_label(sD,’replace’, [51:100],’Ve’); 
sD=som_label(sD,'replace’,[101:150],’Vi’); 
sM=som_make(sD) 

sM=som_autolabel(sM,sD,’vote’); 

colormap(l-gray) ; 7. bedzie w~odcieniach szarości 
som_show(sM,’umat’,’all’,’comp’,1); 

som_show_add(’label’,sM,’TextSize’,8,’TextColor’, ’r’,... 

’subplot’,2); 

VI - 

// job_som2.m 

7. Przeczytanie danych, ich normalizacja 

sD=som_read_data(’iris.data’); sD=som_normalize(sD,’var’) 

"/, Utworzenie struktury—mapy, nadanie etykiet przez ’freq’ 

"/, i~pokazanie dwóch map na razie bez etykiet 
sD = som_label(sD, 'replace’, [1:50]’, ’Se’); 
sD = som_label(sD, 'replace’, [51:100]’, ’Ve’); 
sD = som_label(sD, 'replace’, [101:150]’, ’Vi’); 
sM=som_make(sD) 

sM=som_autolabel(sM,sD,’freq’); l ’add’, ’addl’, ’freq’ 
som_show(sM,'empty’,'Labels 3 Irises’, 'empty’, 'For hits’); 

"/. Pokazanie etykiet na mapie 

som_show_add(’label’,sM,’TextSize’,10,’TextColor’, ’r’,’subplot’,1); 

"/, Obliczenie tzw. trafień ("hits"): 

"/, ile próbek reprezentują poszczególne neurony 

"/o liczba trafień zostanie pokazana na mapie przez 

"/, wielkość heksagonow, a~gatunek irisa odmiennym kolorem 

hl=som_hits(sM,sD.data(l:50,:)); h2=som_hits(sM,sD.data(51:100,:)); 

h3=som_hits(sM,sD.data(101:150, : )); 

som_show_add('hit’,[hl, h2, h3],'MarkerColor’,... 

[100; 010; 00 1], 'subplot’, 2); 

"/, Otrzymujemy plansze z~dwoma mapami: z~etykietami i~trafieniami 

Ti - 
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Rysunek 9.2: Dwa przykłady map Kohonena otrzymane za pomocą pakietu somtool- 
box2. Góra: Mapa u-mat - po lewej, mapa cechy nr 1 (SepalL) po prawej. Dół: Nazwy 
wektorów kodowych po lewej (opcja ’freq), liczby trafień (hits) znakowane wielkością 
heksagonu - po prawej, 


9.5. Wizualizacja mapy — procedura som_grid 

Procedura pełni funkcję podobną do mesh w MatLabie. Przedstawia graficznie dane 
dwu- i trójwymiarowe zawarte w strukturach sD i sM, jednocześnie pozwala na swo¬ 
bodne operowanie markerami, kolorami i grubościami linii dostępnymi w MatLabie. 
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9.5.1. Przykładowy skrypt wykorzystujący procedurę som_gnd 

Pokazujemy użycie procedury sormgrid - instrukcje wypisane z modułu demonstra¬ 
cyjnego som_demo2.m 

7.7.- 

7.7. job_grid.m — som_demo2.m 
sD=som_read_data(’iris.data’); 

7. teraz skrócenie etykiet próbek 
sD = som_label(sD,'replace’, [1:50]’,’se’); 
sD = som_label(sD,'replace’, [51:100]’,’ve’); 
sD = som_label(sD,'replace’,[101:150]’,’vi’); 

7. wykres typu scatterplot-matrix 

7. Here are the histograms and scatter plots of the four variables. 
fhl=figure; k=l; 
for i=l:4, 
for j=l:4, 
if i==j, 

subplot (4,4,k) ; hist (sD. data(: , i) ) ; title (sD. comp_names{i}-) 
elseif i<j , 

subplot(4,4,k); 

plot(sD.data(:,i),sD.data(:,j), ’k.’) 
xlabel(sD.comp_names{i}) 
ylabel(sD.comp_names{j}) 
end 
k=k+l; 
end 

end; drawnow; 

7.7. Normalizacja danych i"utworzenie mapy 

7.7. 

sD = som_normalize(sD,’var’); 
sM = som_make(sD); 

fh2=figure; subplot (2,2,1) ; 7.The map grid in the output space . 
som_grid(sM,'Linecolor’,’k’); view(0,-90), title(’Map grid’) 

xlabel(’m2’); ylabel(’ml’); 

7. - A surface plot of distance matrix: both color and 

7. z-coordinate indicate average distance to neighboring 

7. map units. This is closely related to the U-matrix. 

subplot(2,2,2) 

Co=som_unit_coords(sM); U=som_umat(sM); U=U(1:2:size(U,1),1:2:size(U,2)); 
som_grid(sM,’Coord’,[Co, U(:)],’Surf’,U(:),'Marker’,'nonę’); 
view(-80,45), axis tight, title('Distance matrix’) 
colorbar; xlabel(’m2’); zlabel(’dist’); ylabel(’ml’); 

7. - The map grid in the output space. Three first components 

7. determine the 3D-coordinates of the map unit, and the size 
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"/, of the marker is determined by the fourth component. 

"/, Notę that the values have been denormalized. 

subplot(2,2,3) 

M = som_denormalize(sM.codebook,sM); 

som_grid(sM,’Coord’,M(:,1:3),'MarkerSize’,M(:,4)*2) 

view(-80,45), title('Prototypes ’) 

xlabel('SepalL’); ylabel('SepalW’); zlabel('PetalL’); 

"/, - Map grid as above, but the original data has been plotted 

"/, also: coordinates show the values of three first components 

“/, and color indicates the species of each sample. Fourth 

"/, component is not shown. 

subplot(2,2,4) 

som_grid(sM,’Coord’,M(:,1:3),'MarkerSize’,M(:,4)*2), hołd on 

D = som_denormalize(sD.data,sD); 

plot3(D(l:50,1),D(1:50,2),0(1:50,3),’r.. 

0(51:100,1),D(51:100,2),D(51:100,3),’g.. 

0(101:150,1),D(101:150,2),D(101:150,3),>b.’) 
view(-72,64), axis tight, title('Prototypes and data’) 
xlabel('SepalL’); ylabel('SepalW’); zlabel('PetalL’); 

U 
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SepalL 





PefelW 10 

Pptall 



Rysunek 9.3: Wykres typu scatterplot-matrix dla 4 cech irysa, otrzymany za pomocą 
skryptu job^grid.m jako figurę fhl 




Rysunek 9.4: Fragmenty rysunków wykonywanych przez moduł som-demo2 pakietu 
somtoolbox. Siatka mapy, odległości między węzłami mapy, oraz wykresy trójwymia¬ 
rowe wektorów kodowych oraz danych. Wykonane za pomocą skryptu som_grid 






























Rozdział 10 


Wizualizacja 

wielowymiarowych danych: 
PC A i Odwzorowanie 
Sammona 


Przedstawiamy tu dwie metody wizualizacji punktów znajdujących się w wielowy¬ 
miarowej przestrzeni. Są to: metoda składowych głównych i metoda odwzorowań 
Sammona. Obydwie metody pozwalają dokonać specyficznej redukcji wymiarowości 
do przestrzeni o niższej liczbie wymiarów, w szczególności do d = 2 lub d = 3 wy¬ 
miarów. Przetransformowane dane mogą być wtedy zobrazowane na płaszczyźnie lub 
w przestrzeni 3-wymiarowej przy użyciu zwykłych sposobów. 

Metoda PCA wykonuje transformacje liniowe, a metoda Sammona — transfor¬ 
macje nieliniowe danych. 

Obydwie metody są zrealizowane w pakiecie somtoolbox2. 


10.1. Metoda PCA czyli składowych głównych 

10.1.1. Ogólne zasady wyznaczania składowych głównych 

Metoda składowych głównych jest jedną z bardziej rozpowszechnionych metod ana¬ 
lizy danych wielowymiarowych, szczególnie w kontekście redukcji wymiarowości i wi¬ 
zualizacji. Dobry i zrozumiały opis tej metody znajduje się w książkach Morrisona 
[6] oraz Johnsona i Wicherna [17]; istnieje również świetna monografia lana Jolliffe’a 

W- 
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W środowisku statystycznym metoda składowych głównych jest wiązana z nazwi¬ 
skiem Hotellinga; w środowisku technicznym (przetwarzanie sygnałów, sieci neuro¬ 
nowe) metoda ta jest określana mianem transformaty Karhunen-Loevego. 

Przypuśćmy, że mamy daną tablicę X o wymiarach [n x d\. Wektory-wiersze tej 
tablicy mogą być interpretowane jako punkty £ R d . Przyjmijmy — dla wprowadza¬ 
jącej ilustracji — że d = 2. Wtedy wektory-wiersze tablicy X można zobrazować na 
płaszczyźnie <xl,x2>. Nanosząc dane n punktów (określonych w tablicy X) na tę 
płaszczyznę otrzymujemy tzw. chmurę punktów indywidualnych. Na ogół, dla danych 
statystycznych spotykanych w praktyce, tzn. nie wygenerowanych za pomocą spe¬ 
cjalnych modeli, chmura punktów indywidualnych przypomina elipsę (w przypadku 
d = 3 elipsoidę, a w przypadku d > 3 hiper-elipsoidę). 

Kierunki osi głównych tej elipsy lub (hiper) elipsoidy wyznacza się z równania 
macierzowego 

(Sx - AI)a = 0, (10.1) 

gdzie a = [oi,..., ad] T jest jednym z kierunków głównych, a S.\- jest macierzą kowa¬ 
riancji danych zapisanych w tablicy X. 

Dla jednoznaczności rozwiązania równania (10.1) przyjmiemy, że wektor a jest 
wektorem o długości 1 , co oznacza że a T a = 1 . 

Wektory a spełniające równanie (10.1) i wyznaczające kierunki główne są poszu¬ 
kiwane sekwencyjnie. 

Najpierw szuka się wektora a x wyznaczającego najdłuższą oś elipsy lub (hiper) 
elipsoidy. Wektor ai musi spełniać równanie macierzowe (10.1). Okazuje się, że rów¬ 
nanie to ma d rozwiązań, odpowiadających d wartościom własnym macierzy S. Tak 
więc, rozwiązując równanie (10.1) otrzymujemy d wartości własnych. Porządkując 
otrzymane wartości własne w ciąg nierosnący otrzymujemy 


Ai > A 2 > ... > A d > 0. 


Z każdą z tych wartości własnych jest związany odpowiadający jej wektor własny 

a l, a 2 , ■ ■ ■ , a d- 

Można wykazać, że kierunek najdłuższej osi elipsoidy jest identyczny z kierunkiem 
wyznaczanym przez wektor ai. Dalsze osie elipsoidy są wyznaczane przez dalsze 
wektory własne. 

Wszystkie wektory własne możemy zestawić w macierz A: 

A = [ a !, a 2 ,... ,a d ]. 

Można wykazać, że macierz A wyznacza transformację będącą rotacją układu 
współrzędnych <xl,..., xd> (zaczepionego w punkcie (xl,..., xd)) do układu wek¬ 
torów <ei,... ,e d > wyznaczających osie główne elipsoidy. 
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Rotacja ta może być użyta do transformacji macierzy X. Otrzymamy wtedy: 


Y = XA. (10.2) 

Otrzymana w ten sposób tablica Y przedstawia te same dane co tablica X, ale 
wyrażone w innym układzie współrzędnych. 

Definicja. Elementy tablicy Y określonej wzorem (10.2) są nazywane współrzęd¬ 
nymi głównymi (ang. principal coordinates) tablicy danych X. 

Macierz kowariancyjna „nowych” danych Y równa się 


Sy = 


Ai 

0 


0 

A2 


0 

o 


0 0 ■ • ■ A p 


(10.3) 


Ze wzoru (10.3) wynika, że przetransformowane dane mają następujące własności: 


• wariancje wartości zapamiętanych w kolejnych kolumnach wynoszą odpowied¬ 
nio Ai, A 2 ,..., A p , 

• pary kolumn tablicy Y przedstawiają wartości zmiennych nieskorelowanych. 


Można pokazać, że 


d d d 

trace(S x ) = ^ Jki^ = ^ sf = ^ \j. 

1—1 elementy przek. 1—1 wariancje xl,...,xd 1—1 

Wynika stąd, że suma wariancji współrzędnych głównych równa się sumie wariancji 
zmiennych wyjściowych danych w układzie < xl ,..., xp >. 

Jeżeli suma dwóch lub trzech pierwszych wartości własnych jest duża (np. stanowi 
90% lub więcej całkowitej ich sumy), to znaczy że wariancje wyjściowych zmiennych 
A'i,..., X p są odtwarzane w znacznym stopniu przez nowe (przetransformowane) 
zmienne, czyli przez współrzędne główne. Można stąd dedukować, że elipsoida w R d 
obejmująca dane X jest faktycznie mocno spłaszczona, wobec czego jest możliwe 
i uzasadnione posługiwanie się reprezentacją tych danych w niższej przestrzeni R r 
wyznaczonej przez r pierwszych składowych głównych. 

Zauważmy również, że każda współrzędna główna (składowa główna) jest otrzy¬ 
mywana jako kombinacja liniowa wszystkich zmiennych wyjściowych: 


ykd = X x aj. 


(10.4) 
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10.1.2. Wskazówki praktyczne i podsumowanie 

Przy obliczaniu współrzędnych głównych (składowych głównych) należy obejrzeć ko¬ 
niecznie wartości własne odpowiadające tym składowym. 

Tylko wtedy, gdy udział pierwszych wartości własnych w stosunku do sumy 
wszystkich wartości własnych jest duży, mamy teoretyczne podstawy (i uzasadnienie, 
jesteśmy uprawnieni) do zastąpienia tablicy danych Xod kolumnach zredukowaną 
macierzą Y( r ) o r kolumnach (r < d). 

Podsumowanie 

Współrzędne główne (składowe główne) są to współrzędne w układzie osi głównych 
elipsy lub (hiper) elipsoidy koncentracji przybliżającej chmurę punktów indywidual¬ 
nych w R d . 

Współrzędne te są dane wzorem 

Y = XA , 

gdzie A zawiera kolumnami wektory własne macierzy kowariancji Sx, a X jest scen- 
trowaną tablicą danych. 

Jeżeli — wskutek współzależności analizowanych cech — elipsoida koncentracji ma 
kształt mocno spłaszczony, wtedy jest uzasadnione zastąpienie d współrzędnych ory¬ 
ginalnych (zapisanych w tablicy X) mniejszą liczbą r współrzędnych wziętych jako 
pierwsze r kolumny z tablicy Y = XA. 

Utrata informacji jest mierzona ilorazem 

d d 

E vE A - 

j=r +1 j = 1 

przy czym 

Ai + ••• Ad = $! + ••• + Sd, 

gdzie s 2 j (j = 1, ..., d) oznacza wariancję zmiennej Xj. 


10.1.3. Obliczenia PCA za pomocą pakietu somtoolbox 

Obliczenia składowych głównych odbywają się za pomocą funkcji pcaproj. Przykła¬ 
dowe wywołania tej funkcji: 

(•) [P,V,me,l] = pcaproj(D, odim); gdy liczymy po raz pierwszy, 
lub 

(•) P = pcaproj(D, V, me); gdy znamy już operatory rzutowania. 

Jako argument funkcji pcaproj należy podać przede wszystkim tablicę danych 
którą chcemy zrzutować do układu współrzędnych głównych. Może to być tabli¬ 
ca D według standardu MatLaba; można tu również posłużyć się odpowiednimi 
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strukturami data_struct SD lub map_struct sM — specyficznymi dla pakietu som- 
toolbox. 

(•) Pierwsze, standardowe wywołanie, dostarcza dla zadeklarowanej tablicy D odim 
współrzędnych głównych. Obliczenia procedury są w tym przypadku następujące: 

Najpierw - dla zadeklarowanej tablicy D — oblicza się wektor średnich me oraz 
macierz kowariancji Sx- 

Następnie wyznacza się wszystkie wartości własne i wektory własne macierzy S y- 
Obliczone wartości własne — wraz z odpowiadającymi im wektorami własnymi 
— są porządkowane malejąco. 

Wreszcie oblicza się — według wzoru (10.4) — żądanych odim. współrzędnych 
głównych, które umieszcza się jako wynik w tablicy P o wymiarach [n x odim]. 

Oprócz tego w wynikach mogą być umieszczone: 

o V o wymiarach [d x odim ] - wektory własne użyte do projekcji (pierwszych 
odim kolumn m. A używanej w naszym opisie metody), 

❖ me o wymiarach [1 x d] - wektor średnich, 

o 1 o wymiarach [1 x odim] - frakcje |A j\/J2j |Ay|. 

(•) Obliczone przy pierwszym wywołaniu wektory własne mogą służyć, wraz z wek¬ 
torem średnich, do projekcji innego zbioru danych, np. zapamiętanego w tablicy 
Dl - oczywiście o tej samej liczbie zmiennych, co tablica D z której wyznaczono 
wektory V. 

Rzuty dla innego zbioru danych, np. Dl możemy otrzymać za pomocą rozkazu: 

PI = pcaproj(Dl, V, me) 

10.1.4. Przykładowy skrypt do obliczeń współrzędnych głów¬ 
nych za pomocą pakietu somtoolbox 

Podany niżej skrypt oblicza najpierw 3 pierwsze składowe główne dla znormalizo¬ 
wanych danych zapamiętanych w strukturze sD. Ponieważ dane są znormalizowane, 
obliczona macierz kowariancji jest faktycznie macierzą korelacji rozważanych zmien¬ 
nych. Frakcje odtworzonej sumy X^=i A-i przez kolejne składowe wynoszą odpowied¬ 
nio: 0.7272, 0.2303, 0,0368, a skumulowane sumy tych frakcji: 0.7272, 0.9575, 0.9943 
Najpierw rzutujemy do wyznaczonego układu kierunków głównych wektory ko¬ 
dowe zapamiętane w strukturze irysowej sM. Rzutowane wektory kodowe są znorma¬ 
lizowane. 

Wykres współrzędnych głównych jest robiony za pomocą procedury som_grid 
taki wykres pokazuje powiązania wektorów kodowych odpowiadające sąsiedztwu na 
mapie sM. 
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Na ten sam trójwymiarowy wykres zostają zrzutowane — za pomocą obliczonego 
już operatora rzutowania czyli tablicy V wektory kodowe zapamiętane w struktu¬ 
rze sD. Jest to 150 punktów £ R d odpowiadających poszczególnym kwiatom irysa. 
Gatunki poszczególnych kwiatów są zaznaczone różnymi kolorami: czerwony oznacza 
irysy Setosa (Se), zielony - Versicolor (Ve), a niebieski - Virginica (Vi). Som-quality: 
q=0.3927, t=0.0133. 

Odpowiedni wykres jest pokazany jako Rysunek 10.1. 

Wykres samych wektorów kodowych (znormalizowanych) jest pokazany na ry¬ 
sunku 10.2. Rysunek ten jest sporządzony w tym samym układzie współrzędnych, co 
rysunek 10.1, tj. w układzie składowych głównych pcal,pca2,pca3 wyznaczonych ze 
znormalizowanych danych obejmujących 150 kwiatów irysa zapamiętanych w struk¬ 
turze sD. Som,-qiiality: q=0.3927, t=0.0133. 



Rysunek 10.1: Wykres danych ze struktury sD (punkty zaznaczone trzema kolorami 
w zależności od gatunku irysa) oraz wektorów kodowych ze struktury sM (punkty 
połączone siatką) w układzie składowych głównych pcal,pca2,pca3. Układ współrzęd¬ 
nych obliczony ze znormalizowanych danych irysa 


A oto skrypt, za pomocą którego wykonano oba rysunki: 

"/."/•-job_pca.m- 

°/, make the data and the S0M 

sD = som_read_data(’iris.data’); sD = som_normalize(sD,’var’); 
sD = som_label(sD, ’replace’, [1:50], ’Se’); 
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Rysunek 10.2: To samo, co na rysunku 10.1, ale pokazano tylko siatkę wektorów 
kodowych wraz z ich etykietami 


sD = som_label(sD, 'replace’, [51:100], ’Ve’); 
sD = som_label(sD, 'replace’, [101:150], ’Vi’); 
sM = som_make(sD); sM = som_autolabel(sM,sD,’vote’); 

7. pause 

‘/, find PCA-projection of the data — from correlation matrix 
[Pd,V,me,lambda] = pcaproj(sD,3); 

‘I, coordinates, eigen-vectors, center, lambda 
’/, plot the map grid projection 

7. from help: [P,V,me,l] = pcaproj (D, odim) 

7, or P = pcaproj (D, V, me) 

fhl=figure 

som_grid(sM,’Coord’,pcaproj(sM,V,me),'marker’,'nonę’, ... 

'LineColor’,’c’,'Label’,sM.labels,'LabelSize’,10,’labelcolor’,’k’); 
xlabel(’pcl’), ylabel(’pc2’), zlabel(’pc3’) 

‘/, plot also the original data with color indicating subspecies 
hołd on, grid on 

colD = [repmat([l 0 0],50,1); repmat([0 1 0],50,1);... 
repmat ( [0 0 1] , 50,1) ] ; 

som_grid(’rect’,[150 1], ’Line’, ’none’,’Coord’,Pd,'markercolor’,colD) 

7.7. 

7. make another map containing projections of codebook vectors only 
7. but using the projections obtained from the normalized data set 
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fh2=figure 

som_grid(sM,’Coord’,pcaproj(sM,V,me),'marker’,’none’, ... 

'LineColor’,’c’,'Label’,sM.labels,'LabelSize’,10,’labelcolor’,’k’); 
xlabelC’pcl’), ylabel(’pc2’), zlabel(’pc3’) 
axis([-3, 4, -4, 4, -1.5, 1]); grid on 

i 

u - 


10.2. Odwzorowanie Sammona 

10.2.1. Omówienie algorytmu odwzorowania Sammona 

Mamy n wektorów d-wymiarowych x^ (?’ = 1,2, ...,n) które można interpretować 
jako punkty osobnicze leżące w przestrzeni R d . Chcemy zrzutować te punkty do 
przestrzeni M-wymiarowej (M = 2,3); odpowiednie rzuty będziemy oznaczać jako 
y.i = [j/i,..., um] t ■ Tak więc dla każdego punktu x,; należy znaleźć jego rzut yy 

Xi => yi , gdzie x, ; £ R d , y { £ R M . 

W obu przestrzeniach (tj. R d i R M ) określamy pojęcie odległości d(i,j ) = djj mię¬ 
dzy punktami o numerze i oraz j leżącymi w tej przestrzeni. Do określenia odległo¬ 
ści można zastosować dowolną metrykę, w szczególności euklidesową. W przypadku 
przyjęcia tej ostatniej odległość między dwoma punktami dehniuje się jako pierwia¬ 
stek kwadratowy z sumy kwadratów różnic liczonych po wszystkich składowych tych 
punktów. 

Wprowadźmy następujące oznaczenia: 

D*j = Z?(x,;,Xj) - odległości między punktami w przestrzeni R d , 

= D(yi,yj) - odległości między odpowiednimi rzutami w przestrzeni R M . 

Zadanie odwzorowania nieliniowego Sammona polega na takim doborze wektorów 
y i; aby zminimalizować funkcję błędu E zdefiniowaną w postaci 

1 n 

E =-Y,\ D h- D ijf/ D iv ( 10 - 5 ) 

i<j 


= Y, D b 


i<j 


Dij — 


M 


\ \ El (dis Ujs) ) D tJ — 


-U 'y 1 ( x is x js ) • 


przy czym 
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W minimalizacji funkcji błędu E określonej wzorem (10.5) Sammon zastosował ite- 
racyjną metodę optymalizacyjną Newtona, uproszczoną do postaci 


y pq {k + 1) = y pq {k) - rjA pq (k), 


( 10 . 6 ) 


w której 


A pq (k) — 


SE 
Sy pq 
I S 2 E 

1 s 2 vŁ 


k oznacza numer iteracji, 

A P q(k) jest ilorazem odpowiedniej składowej gradientu i diagonalnego składnika he- 
sjanu — wyznaczonych w k-tej iteracji, 

r] jest współczynnikiem uczenia; najczęściej nadaje się mu wartości z przedziału 
[0.3; 0.4], 

Przy definicji funkcji błędu w postaci (10.5) odpowiednie składowe gradientu 
i hesjanu są dane wzorami: (źródło: Osowski, str. 267 268) 


SE 

Sy P q 


E 


n*. _ d 

U p3 EJpj 


C j=hi*P DpjD Pi 


][Vpq Vjq] j 


(10.7) 


Ó 2 E 

Pyf q 


-i e 


c j=h ^ P D pj D pj 


\(d;, - d„) - 


n* _ n 

u p:i u p 


D 


pj 


D 


')]• ( 10 - 8 ) 


pj 


Odwzorowanie Sammona jest odwzorowaniem nieliniowym punktów z przestrzeni 
R d na odpowiednie ich „rzuty” leżące w przestrzeni R AI . 


10.2.2. Obliczenia odwzorowania Sammona za pomocą pakietu 
somtoolbox i przykładowy skrypt 

Odwzorowaniu mogą podlegać zarówno wektory danych jak i wektory wagowe repre¬ 
zentujące klastery danych. 

Procedura sammon zaimplementowana w pakiecie somtoolbox2 może być stosowa¬ 
na zarówno do tablicy danych D, struktury danych czyli data_struct sD jak również 
do wektorów kodowych znajdujących się w strukturze som.struct sM. 

Przykładowe wywołania procedury: 

P = sammon(D,2); P = sammon(sM,3); P = sammon(sM,3,[],[],[],Md); 

Po obliczeniu tablicy projekcji P może nastąpić wizualizacja graficzna przetrans- 
formowanych danych znajdujących się w tablicy P za pomocą procedury som_grid lub 
też innych procedur wizualizacyjnych MatLaba. 

Przykładowe wywołanie procedury som_grid: 
som_grid(sM,’Coord',P); lub też, z bardziej wyspecyfikowanymi parametrami: 
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Rozdział 10. Wizualizacja wielowymiarowych danych 
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Rysunek 10.3: Znormalizowane Wektory kodowe danych ’Iris’ zrzutowane metodą 
Sammona na dwu-wymiarową płaszczyznę < yl,y2> (M = 2) 


som_grid(sM, , Coord',sammon(sM,2),’LineColor','r , ,'Laber,sM.labels, ... 

' la bel Color ’, ’ b’, 'Marker','d', MarkerColor, ' k') 

Otrzymamy wtedy rysunek 10.3. 

Na rysunku tym obserwujemy zjawisko skręcania się siatki (mapy). Obraz taki 
jest dość częstym zjawiskiem przy stosowaniu odwzorowań Sammona. Wynika on 
z trudności otrzymania obrazu dwuwymiarowego który ma pokazać odległości między 
punktami w przestrzeni wielowymiarowej, a więc w przestrzeni o nieporównywalnie 
bogatszych możliwościach. 

Czasami pomaga wystartowanie z innej mapy (tj. uzyskanej z innego przybliżenia 
początkowego przy uczeniu mapy, lub innego układu neuronów na mapie). 

A oto przykładowy skrypt: 

"/."/•-job_samm.m- 

"/, make the data, shorten labels, make the S0M 

sD = som_read_data(’iris.data’); sD = som_normalize(sD,’var’); 

sD = som_label(sD, ’replace’, [1:50], ’Se’); 

sD = som_label(sD, 'replace’, [51:100], ’Ve’); 

sD = som_label(sD, 'replace’, [101:150], ’Vi’); 

sM = som_make(sD); sM = som_autolabel(sM,sD,’vote’); 

"/,sM = som_make(sD) ; sM = som_autolabel(sM, sD, ’ vote ’ ) ; 

"/, visualization of codebook vectors by Sammon projection 
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’/, Ps=sammon(sM,3) ; 

som_grid(sM,'Coord’,sammon(sM,3),'marker’,'nonę',... 

'Label’,sM.labels,’labelcolor’,’k’); 

*/.- 





Rozdział 11 


Backpropagation — 
propagacja wsteczna 

11.1. Oznaczenia wstępne 

Rozpatrzmy sieć dwuwarstwową, składającą się z warstwy wejściowej (warstwa „0”), 
warstwy pośredniej (ukrytej) „1” i warstwy wyjściowej „2”. 

Niech symbole i, j, k oznaczają numery jednostek (neuronów) należących do tych 
warstw, jak to pokazano na rysunku 11.1. 



aj, Zj = g(aj) 


Rysunek 11.1: Przekazywanie sygnałów od neuronu i na wejściu do neuronu j warstwy 
pośredniej, a stamtąd do neuronu k na wyjście 


Przy określaniu par < i, j > oraz < j,k > obowiązują następujące zasady: 
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Rozdział 11. Backpropagation — propagacja wsteczna 

o Porządek występowania elementów w parach jest ściśle określony. 

o Pierwszy element należy do warstwy niższej (wcześniejszej), drugi do warstwy 
następnej, wyższej. 

o Połączenie między parami jest skierowane; bodźce (sygnały) są przekazywa¬ 
ne od jednostki stojącej na początku pary (tj. jednostki z niższej warstwy) 
do jednostki stojącej na końcu pary (tj. jednostki umiejscowionej w wyższej 
warstwie). 

Wymienione fakty można prześledzić na rysunku 11.1. 


11.2. Przekazywanie bodźców 


Popatrzmy teraz na elementy i,j , k pod kątem ich aktywacji i przekazywanych dalej 
bodźców. 

Neuron i należący do warstwy „0” ma aktywację = X{. 


Neuron j zbiera bodźce od I jednostek z warstwy poprzedniej i sumuje je z odpo¬ 
wiednimi wagami. Oznacza to obliczanie sumy 


dl) 


2 >« 4 0) . 


gdzie uj^I 1 oznacza wagę połączenia skierowanego od jednostki i do jednost¬ 
ki j; wskaźnik d) u góry oznacza, że są to wielkości obliczone przez neuron 
umiejscowiony w warstwie „1”. 


Na podstawie tak wyznaczonej sumy zostaje wyznaczona aktywacja j-tego 
neuronu 


z j 1} = 


gdzie g(.) oznacza funkcję aktywacji przyjętą dla warstwy „1”. 

Obliczona aktywacja stanowi bodziec (sygnał) przekazywany dalej (for- 
wards, czyli w przód) do następnej warstwy, w szczególności do k- tej jednostki 
tej warstwy (patrz rysunek 11.1). 


Neuron k zbiera bodźce od J jednostek pochodzących z warstwy poprzedniej (na 
rysunku 11.1 ma ona numer „1”) i sumuje je z odpowiednimi wagami. Oznacza 
to obliczenie sumy 


4 2) = £«£>•#’, 


fce/c 


gdzie w^j 1 oznacza wagę połączenia w parze < j, k >, a wskaźnik ^ 2 ) u góry 
oznacza, że są to wielkości obliczone przez neuron umiejscowiony w warstwie 
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(2) 

Na podstawie tak wyznaczonej sumy aj“ zostaje wyznaczona aktywacja fc-tego 
neuronu 


( 2 ) _ 
k 



), 


gdzie g{aj, J ) oznacza funkcję aktywacji przyjętą dla warstwy „2”. Może to być 
ta sama funkcja, którą przyjęto dla warstwy „1”, w takim przypadku g(-) = g(-). 

Jeżeli neuron nr k jest neuronem warstwy wyjściowej, to obliczona aktywacja 
z) jest wynikiem końcowym y^ i jest przekazywana na wyjście jako wynik y 


Rezultaty tych działań są podsumowne w tabelce poniżej: 


Warstwa Neurony 


warstwa „2” neuron O k 

\T 

warstwa „1” neuron O j 

\T 

warstwa „0” neuron Q ' 


Formuły obliczeniowe 


( 2 ) 

a k =Hk W kjZ; , 


( 1 ) 


4 2 ) = 5 ( 4 2) )> Vk=Zk ] 


(1) t 

a ) = Ej WjiĄ 


(0) 




(o) 

z\ = Xi 


Tabela 11.1: 


11.3. Określenie wektora wag 


W omawianym wyżej modelu wagi występują w 2 warstwach: W warstwie pierwszej 
mamy wagi wynikające z połączeń I jednostek warstwy poprzedniej z j-tą jednostką 
(j = 1,... ,J) obecnie analizowanej warstwy pierwszej 


J 1 ) 


r (!) 


, • (1) r C 1 ) (l)lT 

'i > • ■ • > ™ j- , gdzie w y = [w) 1 ',...,w) I '\ . 

Analogicznie, dla warstwy drugiej mamy wagi 


W (2) W (2) 

W 1 1 • • • 1 w /i' I 


gdzie wj . 2 ' 1 = [u; 


,,( 2 ) 


„(2)lT 


L, fcl ) • • • > ^kj J 


Wszystkie te wektory wag możemy zestawić w jeden wspólny wektor w określony 
następująco: 


= [(wS 1, ) T ,., (w™) 7 ', (w< 2| ) T ,.... (wg>f] 
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Rozdział 11. Backpropagation — propagacja wsteczna 


11.4. Funkcja błędu 


Funkcja błędu E wyznaczona dla N próbek jest zdefiniowana jako suma błędów 
zaobserwowanych dla każdego wzorca n, (n = 1,, N): 


E = ^ E n 


n 


gdzie E n jest błędem dla wzorca n. 

Błąd E n może być definiowany na rozmaite sposoby. Jedną z częściej stosowanych 
definicji błędu jest definicja oparta na zasadzie najmniejszych kwadratów, tj. 



gdzie = y£(x n ,w) jest wynikiem wyprodukowanym przez sieć (sygnałem post- 
synaptycznym) na wyjściu k jako odpowiedź na dany sygnał x" (sygnał presynap- 
tyczny) - przy użyciu wektora wag w charakteryzujących siłę połączeń między neu¬ 
ronami. 

Błąd E n dla n-tego wzorca możemy uszczegółowić, wyrażając go za pomocą pa¬ 
rametrów sieci przedstawionych na rysunku 11.1: 



E n 



E, 


Jak wiadomo, proces uczenia przebiega w ten sposób, że na podstawie prezentowa¬ 
nych wzorców sieć adaptuje (koryguje) swoje wagi tak, aby uczynić ostateczny błąd 
E możliwie najmniejszym. 

W jaki sposób sieć ma to zrobić, tzn. jak skorygować swoje wagi aby osiągnąć 
minimum funkcji błędu? 

Na ogół, jeśli funkcja jest różniczkowalna, robi się to metodami gradientowymi. 
Metody te wymagają obliczenia pochodnych minimizowanej funkcji. 
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11.5. Obliczanie pochodnych funkcji błędu względem 
wag Wji oraz Wkj 


Wprowadźmy oznaczenia 


= 


dE r 


r(2) 


da 


(i) : 


dE n 

da k ] 


Obliczanie pochodnych względem wag drugiej warstwy 


dE n 

W] 


dE n 



4 2) 


da 


( 2 ) 


dw 


( 2 ) 

kj 



gdzie 


4 2) = 


dE n 


da 


( 2 ) 


l (2) 


ponieważ y k = g(a 


( 2 ) 

fc 


)• 


dE n dy k 
dyk da { k ] 


dE n dg(a^) 
dyk da i 2) 


Obliczanie pochodnych względem wag pierwszej warstwy 

dE n _ dE n da ? _ S (D jo) 


( 11 . 1 ) 


( 11 . 2 ) 


(11.3) 


gdzie ćj 1 ' jest dalej obliczana jako pochodna złożona, uwzględniając fakt, że zmienna 

/o\ /o\ 

a^ ' występuje w każdym składniku a y k ; które miało połączenie z jednostką (neuro¬ 
nem) j warstwy pośredniej. 

Aby obliczyć 5j obliczmy ją jako pochodną złożoną, ponieważ E n zależy od aj za 
pośrednictwem a k . 


s? 1 = 


dE n 


da 


(i) 


„(n 


= £' 


dE n da 


(2) , , 

k ó(. 2) da 


da 


( 2 ) 

k 

( 1 ) 


(11.4) 


Należy teraz obliczyć {da k )/{da^'). 

Przypomnijmy, że a^ = ^2 j Wkjg(af' > )] wobec tego (da^)/(daf' > ) = wff g 1 (aj). 
Tak więc 
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Rozdział 11. Backpropagation — propagacja wsteczna 


(11.5) 

Formuła ta określa propagację wsteczną błędów. 

Aby obliczyć błąd 5j dla jednostki j należy znać błędy 5^ warstwy nadrzędnej. 
Wniosek końcowy: 

Pochodne względem wag są funkcjami odpowiednich delt i aktywacji neuronów. 
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